1. Introduction

Atherosclerosis is a chronic inflammatory disease [37]. It involves several different risk factors [11] that contribute to the whole inflammatory process [10,33]. It starts typically with the penetration of the low density lipoproteins (LDL) into the intima layer of the blood vessel where they are oxidized (ox-LDL). Being a dangerous agent, the ox-LDL triggers an immune response that initiates with the recruitment of monocytes, which are later transformed into macrophages. The macrophages eliminate the ox-LDL and are transformed into foam cells, that along with several other substances form the so-called lipid core. Smooth muscle cells migrate and form a fibrous cap which covers and isolates the lipid core from the blood circulation. The lipid core along with the fibrous cap constitute what is called the atheromatous plaque [35]. This plaque changes the geometry of the blood vessel's lumen by narrowing it and interacts with the blood flow. As the plaque continues to grow, the shear stress of the blood flow through the decreasing cross section of the lumen increases. This force may eventually cause the rupture of the plaque, resulting in the formation of thrombus, and possibly a stroke [20]. This paper is devoted to the FSI analysis of the interactions of the blood flow with the plaque and the vessel wall.

The methods to describe the interaction between fluid and structure are of great importance to achieve accurate numerical simulations. These methods include among others, the level set method [4], the immersed boundary method [31,32], the fictitious domain method [17,18], the Fully Eulerian formulation [14,41] and the Arbitrary Lagrangian Eulerian (ALE) approach [21,9,30]. In this paper we use the ALE method to simulate our FSI model.

A large number of biomechanical and imaging studies started to emerge since the early 2000s addressing the role of the hemodynamic shear stress in the destabilization of vulnerable plaques, investigating the conditions of plaque rupture and using a FSI problem between blood flow and atherosclerotic plaque [27,28,40,25]. In the same framework, the coupling between the structure and the fluid was applied to FSI problems in which the coupling was evaluated at the blood flow and the vessel wall interface. In order to capture the pulse wave propagation due to the interaction between blood and the vessel wall, a 3D shear-thinning viscosity model to account for blood rheology was proposed in [23], where the blood vessel was represented as an idealized cylinder. Also, in [8], while dealing with a realistic geometric model of the aorta, the authors tested different types of defective BCs to accurately simulate the pressure waves.

In addition to the important choice of a suitable FSI numerical method, as described above, a correct selection of appropriate mathematical models to handle the blood rheology, the mechanical behavior of the vessel wall and the atheromatous plaque is also an important issue. Indeed, the dynamics of blood flow is mainly influenced by its rheological properties and scientists still face the controversy to consider blood as a Newtonian fluid or as a non-Newtonian fluid with inelastic shear-thinning viscosity, and possibly with viscoelastic or yield-stress properties. To investigate the flow effects induced by different blood constitutive equations, four non-Newtonian blood models, namely Power Law, Casson, Carreau-Yasuda, and Generalized Power Law, as well as the Newtonian model were used in [22] to compare the corresponding blood flow behavior and quantify the differences between them to analyze their significance. Another controversy is related to vessel wall and atheromatous plaque modeling. The vascular wall has a very complex nature and deriving an accurate model for its mechanical behavior is rather difficult [34]. Furthermore, arterial walls of major arteries are composed of several layers, each of them with different mechanical characteristics [29]. Moreover, it has been shown that influential geometrical parameters such as the fibrous cap thickness, the stenosis ratio, the lipid core width and length, all combined, determine the plaque vulnerability [6]. Recently, the influence of variability in material properties of media, fibrous cap, lipid and intraplaque hemorrhage/thrombus, obtained from direct material tests on stress and stretch conditions within the plaque structure, was assessed [43]. If the plaque is not calcified, it grows, increasing its volume. A recent mathematical model and simulation combining the plaque evolution with FSI has been presented in [42]. Nevertheless, throughout most of the numerical studies, plaque components were assumed to be incompressible and nonlinear because ideally the human tissue is hyperelastic [27].

After this brief description of the difficulties regarding the modeling of each physical media (blood, arterial wall and atheroma plaque), it is clear that modeling the FSI of the blood-arterial wall and blood-plaque simultaneously is a challenging task.

One can find an attempt in [2], where, after a 2D analytical study to prove the existence of a weak solution for incompressible non-Newtonian fluids with non standard BCs, the authors tested the model numerically by coupling it with a nonlinear hyperelastic model for the arterial wall and the atheromatous plaque, with different values of Young modulus corresponding to each structure. They showed that considering the classical models of interaction between the blood flow (modeled as a Newtonian fluid), the vessel wall and the atherosclerotic plaque, the risk of rupture is less predictive because it underestimates the stress over the plaque comparing with the case where blood is modeled as a non-Newtonian fluid. These previous results were extended by considering a 3D geometry of an idealized blood vessel in order to be more realistic, the aim being to analyze the evolutionary non-linear 3D flow characteristics of blood in a stenosed artery [3]. However, for the numerical simulations based on the theoretical study, the fluid and the structures were described using a 2D axi-symmetric model in order to reduce the overall computational complexity. The parameters describing the blood rheology, the vessel wall, and the plaque were taken from the literature and experimental data. A special attention was paid to the effects of the wall motion on the local fluid displacement, on the stresses, and on strains in the diseased arterial wall.

It is known that the circulatory system can be seen as an infinite extremely complex network including the heart and blood vessels. Therefore, mathematical and computational modeling of blood flow dynamics is a complex problem that requires simplifying assumptions of the physical properties and BCs. Rather than modeling the entire circulatory system, based on full 3D models, which is beyond the capability of current computers, segments of the circulation are studied to reduce the computational time, requiring appropriate inflow and outflow BCs. Several approaches based on data assimilation [19], defective BCs [8] and absorbing conditions [24] have been used by several authors to address this question.

In this study a 3D FSI model incorporating fluid-structure interaction (blood-plaque and blood-wall) as well as non-Newtonian effects in a 3D geometry to analyze the flow and structures stresses is considered. It consists of a generalization of the study done in [3] with more realistic pressure input profile [1] at the inlet and a linear absorbing condition (LAC) at the outlet [24,36].

---

2. Mathematical modeling and numerical methods

---

2.1. Geometry of the atherosclerotic blood vessel

The computational domain is shown in Figure 1. It consists of a cylindrical tube representing the blood vessel with a narrowing of its interior to mimic the stenosis. The flow channel with radius $R=0.3$ cm and height $H=2$ cm contains the fibrous cap and lipid pool described by the trigonometric functions corresponding to a stenosis length of $l=1$ cm. The degree of arterial stenosis is $\approx70\%$ with a fibrous cap thickness of $0.5$ mm. The thickness of the whole wall was chosen to be $1$ mm. To avoid the entrance and exit effects, the computational domain was extended upstream (2 cm) and downstream (4 cm) the stenosis.

---

2.2. Fluid equations

The blood flow is governed by the generalized Newtonian equations for incompressible fluids:

$\begin{array}{c} {\displaystyle\rho_{f}\frac{\partial {\bf{u}}}{\partial t}+\rho_{f}({\bf{u}}\cdot\nabla){\bf{u}}-\nabla\cdot(2\mu({\bf{s}}({\bf{u}})){\bf{D}}{\bf{u}})+\nabla p=0} \\ {\nabla .{\bf{u}} = 0.} \end{array}$

(1)

This evolution problem involves the blood velocity ${\bf{u}}=({\bf{u}}_1, \mbox{ }{\bf{u}}_2, \mbox{ }{\bf{u}}_3)$ and the pressure $p$ as unknowns. The dynamic viscosity $\mu$ is a function of the strain rate tensor ${\bf{D}}{\bf{u}}=\frac{1}{2}(\nabla{\bf{u}}+\nabla^{T}{\bf{u}})$, or more precisely is a function of ${\bf{s}}({\bf{u}})$, the second invariant of the strain rate tensor, defined by:

$({\bf{s}}({\bf{u}}))^{2}=2{\bf{D}}{\bf{u}}:{\bf{D}}{\bf{u}}=2\sum\limits_{i, j}({\bf{D}}{\bf{u}})_{ij}({\bf{D}}{\bf{u}})_{ji}.$

We consider $\mu$ defined by the Carreau-Yasuda law,

$\mu({\bf{s}}({\bf{u}}))=\mu_\infty+(\mu_0-\mu_\infty)(1+(\lambda{\bf{s}}({\bf{u}}))^2)^{(n-1)/2},$

with $\mu_{0}=0.0456$ Pa.s, the viscosity at the lowest shear rate, $\mu_{\infty}=0.0032$ Pa.s, viscosity at the highest shear rate, $\lambda=10.03$ s, $\rho_{f}=1060$ Kg.$\mbox{m}^{-3}$, the fluid density, and $n=0.344$ (corresponding to a shear-thinning viscosity fluid). The values of the parameters are taken from [15]. In the case where blood flow is modeled as a Newtonian fluid, its constant viscosity is given by $\mu=0.0035$ Pa.s. The pressure waveform at the inlet is shown in Figure 2. At the outlet we consider typical natural BCs where we fix the stress

$2\mu({\bf{s}}({\bf{u}})){\bf{D}}{\bf{u}}\cdot{\bf{n}}-p{\bf{n}}=-\bar{p}{\bf{n}}.$

(2)

Although it is known that this BC can affect the velocity profile, these spurious effects are localized on the outlet surface and, therefore, far from our region of interest. As it was shown in [24], these BCs, when used for a FSI simulation, fail to represent the wave propagation phenomena that would be obtained if an extended domain was considered downstream the artificial boundary. In fact, after a single pulse overpasses the artificial boundary, condition (2) introduces high amplitude wave reflections upstream that region (see [24]). To minimize such spurious effect, in [12,24] the authors coupled the 3D model with a 1D model representing the downstream domain. A less expensive method, the so called LAC was suggested in for non-Newtonian FSI blood flow simulations [24], and used successfully in [36] for realistic geometries. Both methods were shown to reduce significantly, yet not totally, the spurious oscillations. Here, we adopt the less expensive method, the use of LAC conditions, given by the expression $\displaystyle\bar{p}=\frac{\sqrt{\rho_f\beta}}{\sqrt{2}A^{5/4}_{0}}Q,$ where $Q$ is the flow rate, $A_0$ the cross-section reference area at rest. The coefficient $\beta$ is related to the mechanical properties of the vessel wall through the expression $\displaystyle\beta=\frac{\sqrt{\pi}h_0E}{1-\nu^2}$, where where $E$ is the Young modulus, $h_0$ is the wall thickness and $\nu$ is the Poisson ratio.

---

2.3. Structure equations

In addition to the complexity of the arterial tissue and the difficulty in obtaining in vivo data, modeling the atheromatous plaque is a heavy task, due to its composition and morphology. In fact, on the one hand, carotid atherosclerotic plaques are composed of lipid deposits with variable numbers of macrophages and T lymphocytes, and the overlying fibrous cap is rich in smooth muscle cells and proteoglycans [39]. On the other hand, plaques with a thick fibrous cap and without necrosis are considered stable, whereas plaques with a thin fibrous cap and areas of necrosis with an increased number of inflammatory cells are considered vulnerable. Also, plaques in which the thrombi occur as a result of complications such as fibrous cap disruption or ulceration are considered unstable [38]. To simplify, the arterial tissue and the fibrous cap are modeled as compliant isotropic materials, ignoring their different composition and fibrous nature. We consider a 3D nonlinear model of hyperelasticity [5] governed by the following equations:

• the equation of motion

$\displaystyle\rho_{s}\frac{\partial^{2}\boldsymbol\eta}{\partial t^{2}}-\nabla\cdot\boldsymbol\sigma ={\bf{F}}, \mbox{ }\boldsymbol\sigma =J^{-1}{\bf{F}}\mathbf{S}{\bf{F}}^T, \mbox{ }{\bf{F}}={\bf{I}}_3+\nabla\boldsymbol\eta, \mbox{ }J=\mbox{det}({\bf{F}}),$

(3)

• strain-displacement equations $\displaystyle\boldsymbol\varepsilon=\frac{1}{2}(\nabla\boldsymbol\eta+(\nabla\boldsymbol\eta)^{T}+(\nabla\boldsymbol\eta)^{T}\nabla\boldsymbol\eta),$ $\boldsymbol\eta$ being the displacement vector, $\boldsymbol\sigma$ the Cauchy stress tensor, ${\bf{F}}$ the deformation gradient tensor, $\mathbf{S}$ the second Piola-Kirchhoff stress tensor, and ${\bf{I}}_3$ is the identity matrix.

The stress-strain relationship is given by

$\displaystyle \mathbf{S}-\mathbf{S}_0={\bf{C}}:(\boldsymbol\varepsilon-\boldsymbol\varepsilon_0-\boldsymbol\varepsilon_{inel}),$

where we denote by $\mathbf{S}_0$ and $\boldsymbol\varepsilon_0$ the initial second Piola-Kirchhoff stress and strain tensors, $\boldsymbol\varepsilon_{inel}$ the inelastic strain tensor, ${\bf{C}}$ the fourth-order stiffness tensor and ":" represents the inner product of two second-order tensors.

In the present study the vessel wall is considered isotropic in such a way that the stiffness tensor has no preferred direction and the corresponding tensor components are computed internally according to:

${\bf{C}}^{ijkl}=\lambda_s\boldsymbol g^{ij}\boldsymbol g^{kl}+\mu_s(\boldsymbol g^{ik}\boldsymbol g^{jl}+\boldsymbol g^{il}\boldsymbol g^{jk}),$

where $\displaystyle\lambda_{s}=\frac{\nu E}{(1-2\nu)(1+\nu)} \mbox{ and } \displaystyle\mu_{s}=\frac{E}{2(1+\nu)}$ are the Lamé coefficients, $E$ is the Young modulus, $\nu$ is the Poisson ratio and $\boldsymbol g$ the metric tensor. Here we consider $\nu=0.45$, $E=3\cdot10^5$ Pa in the upstream and downstream extended sections and $E=6\cdot10^5$ Pa in the vessel wall surrounding the plaque [26], since it has been shown that the presence of a space-occupying the plaque is commonly associated with the expansion of the arterial wall [16]. For the atheromatous plaque, the second Piola-Kirchhoff stress tensor $\mathbf{S}$ is computed as $\displaystyle\mathbf{S}=\frac{\partial\omega}{\partial\boldsymbol\varepsilon}$ where $\omega$ is the elastic strain energy density. According to the Mooney-Rivlin model for hyperelastic materials, choosen to model the lipid pool and the fibrous cap, the energy function reads

$\omega=C_{10}(\bar{I}_1-3)+C_{01}(\bar{I}_2-3)+\frac{1}{2}\kappa(\bar{J}-1)^2,$

(4)

where $\bar{J}$ is the ratio between the current and original volumes and $\bar{I}_1$, $\bar{I}_2$ are the two first strain invariants, modified to be independent of the volume change $(\bar{I}_1=I_1/\bar{J}^{2/3}\mbox{ and }\bar{I}_2=I_2/\bar{J}^{4/3}$, where $I_1$ and $I_2$ are the two first strain invariants). The use of these modified invariants and the last term in equation (4) allows the description of nearly incompressible materials. The parameters $C_{10}$, $C_{01}$ and $\kappa$ are choosen in agreement with experimental measurements, and their values for the fibrous cap and the lipid pool are taken from [28,40] (Table 1):

---

2.4. ALE formulation

Since the vessel wall deforms under the action of the fluid, the mathematical formulation should allow points in the vessel structure to move at every time instant, as well as the points inside the fluid domain. Neither the Lagrangian or the Eulerian formulations allow this behavior without perturbing the artificial upstream and downstream boundaries. Therefore, an ALE formulation is used. We denote by $\boldsymbol w$ the domain velocity. If $\boldsymbol w$ is zero, the method is simply an Eulerian approach, otherwise, if $\boldsymbol w$ is equal to the velocity of the fluid then the approach is Lagrangian. $\boldsymbol w$ can vary arbitrarily and continuously from one value to another in the fluid field.

We rewrite the equations of the fluid, the generalized Navier-Stokes equations, in the ALE formulation as

${\displaystyle\rho_{f}\frac{\partial {\bf{u}}}{\partial t}|_X+\rho_{f}(({\bf{u}}-\boldsymbol w)\cdot\nabla){\bf{u}}-\nabla_x\cdot(2\mu({\bf{s}}({\bf{u}})){\bf{D}}{\bf{u}})+\nabla_xp=0}$

(5)

where $X$ is the Lagrangian coordinate and $x$ is the Eulerian coordinate. From a mechanical point of view, the coupling between the systems (5) and (4) is realized by imposing at the interface:

• the continuity of the velocity over time:

${\bf{u}}=\frac{\partial\boldsymbol\eta}{\partial t}$

• and the equilibrium of the stresses

$-(2\mu({\bf{s}}({\bf{u}})){\bf{D}}{\bf{u}}-p{\bf{I}}_3){\bf{n}}=\mathbf{S}(\boldsymbol\eta)\cdot{\bf{n}},$

where ${\bf{n}}$ is the outward unit vector to the physical interface. This condition should be rewritten on the interface at $t=0$ to be prescribed on the structure model. Using the Piola transform we have

$-(\mbox{det}\nabla_0\boldsymbol\eta)(2\mu({\bf{s}}({\bf{u}})){\bf{D}}{\bf{u}}-p{\bf{I}}_3)(\nabla_0^{-T}\boldsymbol\eta)\cdot{\bf{n}}_0=\mathbf{P}(\boldsymbol\eta)\cdot{\bf{n}}_0,$

where $\mathbf{P}=\mathbf{P}(\boldsymbol\eta)$ is the first Piola-Kirchhoff tensor, $\nabla_0$ indicates the gradient with respect to the Lagrangian coordinates, and ${\bf{n}}_0$ the outward unit vector to the interface at $t=0$ (see [23])

The resulting coupled system of motion and fluid equations was solved using a finite element space discretization based on P1-P1 stabilized elements for the fluid, and P2 elements for the structure. The time discretization was based on the BDF of order 2, with a maximum time step of $5\times10^{-3}$s. The nonlinearities were tackled using Newton's method and each linear iteration solved with the direct PARDISO solver. The simulations were made using Comsol Multiphysics 4.3b [7] in a workstation with twelve Intel® Xeon(R) CPU ES-2620 v2 @ 2.10GHz processors.

In order to improve the accuracy while minimizing the computational cost, we set different element sizes for the mesh approximating the 3D geometry. Therefore, we used smaller element sizes for the lipid pool and the fibrous cap, which are regions characterized by strong gradients and high stresses, than for the vessel wall. Also, the fluid domain was divided into four parts corresponding to different element sizes (see Figure 3), corresponding to a better accuracy in the stenosed region.

All the computations and results presented in the next section were performed with the mesh shown in Figure 3, which corresponds to a total number of degrees of freedom (DOF) of 449 048. Additional simulations for a non-Newtonian Carreau-Yasuda fluid using a coarse and a fine mesh have been performed to analyze the convergence of the numerical solutions with respect to the mesh size (Figure 4). The displacement of the wall was prescribed to be zero only at the artificial boundaries, close to the inlet and outlet fluid boundaries. In the following we call this the moving wall case.

In Figure 5 one can notice that the total displacement curves, in the three meshes, are almost correlated. An exception to this pattern could be the displacement profiles observed after t=0.4s. Although it seems to be an evidence for the oscillations described in section 2.2, which are due to the artificial boundaries, it is not clear why the oscillations increase for the finer mesh. Future analysis, requiring an extra computing power, should include the finer meshes in order to clarify this point. For the velocity magnitude curves, a significant gap is noticed between the coarse mesh curve and the other curves during the time interval corresponding to the peak values of the inlet pressure waveform (see Figure 2). Thus, due to the high computational cost associated to the fine mesh, and that the convergence of the numerical solutions using the coarse mesh is not guaranteed during all the simulation time, an intermediate size mesh was used in this study.

---

3. Results and discussion

---

3.1. Newtonian vs non-Newtonian blood flow

The choice between a Newtonian or a shear-thinning non-Newtonian inelastic model to represent the blood flow is mainly related with the aim of the study. Here we performed numerical simulations for both models in the case of rigid and moving walls to compare the results. Figure 6 gives a comparison between the four cases displaying the total volume displacement of the fibrous cap. One can notice that, practically, there is no significant difference in the total displacement values between the Newtonian model and the non-Newtonian one.

Another important hemodynamic indicator is the wall shear stress (WSS), which is the tangential stress that the fluid exherts on the wall, consisting on the tangential component of the stress tensor, which is given by $\boldsymbol\sigma =p{\bf{I}}-\boldsymbol\tau$, on the wall:

$WSS=\boldsymbol\sigma _n-(\boldsymbol\sigma _n\cdot{\bf{n}}){\bf{n}}=\boldsymbol\tau_n-(\boldsymbol\tau_n\cdot{\bf{n}}){\bf{n}},$

where ${\bf{n}}$ is the outward unit normal to the wall surface, and $\boldsymbol\sigma _n$ and $\boldsymbol\tau_n$ are the normal components of the stress and extra-stress tensors, respectively [15]. The variation of this quantity was observed on the fibrous cap surface. Figure 7 confirms the above observation on the difference between the use of the Newtonian and the non-Newtonian models for blood. Therefore, in the following discussion of the numerical results, we will focus on the 3D FSI model with the fluid described by the Carreau-Yasuda non-Newtonian model. This choice is motivated by the fact that the non-Newtonian behavior is more realistic in its description of the blood flow, since the aim of this paper is to study a 3D FSI model in an idealized stenosed vessel with an absorbing BC imposed directly on the outflow section and its influence on the WSS distribution in the case of moving and fixed walls.

---

3.2. Results for fixed vs moving arterial wall

The numerical results show that the total volume displacement of the fibrous cap at the peak of the pressure ($t$=0.22s) is different between these two cases (0.34 mm and 0.53 mm as maximal displacement values for a fixed and moving walls, respectively). The results also indicate a significant difference for the WSS around the fibrous cap observed in the case of a fixed wall, when compared to the one of a moving wall (Figure 7).

In what follows, we will investigate the distribution of the von Mises stress and WSS over the plaque, since the critical stress indicators are more sensitive to changes of physiological conditions and may provide a better assessment for plaque vulnerability.

Figure 9 plots the von Mises stress in the fibrous cap surface upstream and downstream the stenosis. The stress induced within the plaque was significantly affected by the motion of the wall. In the case of a moving wall, the highest von Mises stress (13.48 KPa) is observed in the region upstream the stenosis, whereas for the fixed wall configuration the highest value (10.45 KPa) is observed downstream the stenosis. This result is in contradiction with the study presented in [3] where the maximum value of the von Mises stress is located upstream the stenosis in the case of the rigid wall. However, two major considerations distinguish our study: (1) the inlet pressure wave and the absorbing BC chosen for the inlet and outlet, respectively and (2) the differentiation between Young modulus values for the vessel wall ($E=3\cdot10^5$ Pa in the upstream and downstream extended sections and $E=6\cdot10^5$ Pa in the vessel wall surrounding the plaque). Furthermore, Figure 8 illustrates how the von Mises stress average changes over the interface of the fibrous cap, in accordance with its distribution over the plaque, as shown in Figure 9. An insignificant difference between the upstream and downstream regions of the stenosis is observed in the case of a fixed wall, but a significant difference exists in the case of moving wall. In terms of interpretation based on the average values, the risk of arterial damage and/or rupture is overestimated for a moving wall, and in the case of a fixed wall, the plaque provides significant load-bearing capabilities.

The WSS distribution on the fibrous cap is depicted in Figure 10. From the plots of different models one can notice that the maximum values upstream (173.27 Pa) and downstream (179.54 Pa) the stenosis in the case of a fixed wall are slightly higher than those observed in the case of a moving wall. However, Figure 11 indicates that the highest average WSS (24.27 Pa) is localized upstream the stenosis with fixed wall, while the difference is insignificant downstream the stenosis between the moving and fixed wall cases. Therefore, the model with rigid wall underestimates the risk of the accumulation of fatty material, since it is commonly believed that intimal thickness will be greater in low shear regions.

---

4. Conclusions

In this work, 3D FSI numerical simulations for blood flow in an atherosclerotic artery have been performed. The computational FSI model considered blood as a shear-thinning viscosity fluid and hyperelastic structures to describe the human carotid artery with an atheromatous plaque. A LAC has been imposed as outlet BC to eliminate the reflexion of pressure waves. To enhance the importance of a fully FSI model, simulations were also performed in the case of a fixed wall structure. The results indicate that the von Mises stress is highly underestimated in the plaque region upstream the stenosis, if a fixed wall model is considered.

Some important aspects still remain to be addressed in order to obtain a more realistic model. Among them we can mention multilayered structure and anisotropic properties of the vessel wall, cell activities, plaque growth and remodeling. Despite these simplifications, the present model demonstrates the importance of the vessel wall movement in terms of the plaque sensitivity to the stress, which may have significant implications for plaque vulnerability assessment.

---

Acknowledgments

FCT (Fundação para a Ciência e a Tecnologia, Portugal) via the grant SFRH/BPD/66638/2009, the project EXCL/MAT-NAN/0114/2012 and the Research Center CEMAT-IST are gratefully acknowledged.

---