Repositorio Institucional de Documentos

Centro Politécnico Superior, Universidad de Zaragoza Trabajo Fin de Máster Máster en Mecánica Aplicada Posgrado en Ingeniería Mecánica y de Materiales Desarrollo de un modelo de daño basado en la microesfera para tejidos biológicos fibrados Autor: D. Pablo Sáez Viñas Director: Dr. Miguel Ángel Martínez Barca Área de Mecánica de Medios Continuos y Teoría de Estructuras. Departamento de Ingeniería Mecánica. Curso 2009/10 Zaragoza, Septiembre de 2010 DESARROLLO DE UN MODELO DE DAÑO BASADO EN LA MICROESFERA PARA TEJIDOS BIOLÓGICOS FIBRADOS RESUMEN Este trabajo fin de máster presenta un modelo de daño micro estructural para tejido vascular usando técnicas multiescala. El trabajo se centra en el comportamiento anisótropo de los tejidos biológicos fibrados y en especial en procesos de ablandamiento o daño. Los tejidos biológicos tienen una composición compleja y muy variable. Algunos de ellos, los denominados tejidos biológicos blandos (vasos sanguíneos, globo ocular, etc.), suelen contener altas concentraciones de colágeno. Esta sustancia se presenta normalmente en forma de fibras, responsables de soportar en gran medida las cargas mecánicas que se producen sobre el tejido. Cuando éstos se deforman lejos de su estado fisiológico se produce el daño del mismo. Con el objeto de estudiar este proceso se ha utilizado una aproximación micro estructural, o más concretamente un modelo basado en la microesfera, para incluir el comportamiento de las fibras. Se ha utilizado un modelo hiperelástico, definiendo una función energía de deformación (FED) en su forma desacoplada. La distribución de las fibras se ha incorporado a través de dos funciones de probabilidad, para tener en cuenta la distribución de las fibrillas alrededor de la orientación preferencial de las mismas. El comportamiento mecánico de cada una de las micro-fibras se ha modelado con dos FEDs diferentes. Además, se ha incorporado un modelo de daño a esta técnica de homogeneización, a través de una formulación termodinámica consistente, que se acopla directamente al modelo hiperelástico de las micro-fibras, para obtener el daño o ablandamiento del tejido. Gracias al modelo micro mecánico, cuando las fibras del material son deformadas, el ablandamiento del material sucede de una manera gradual gracias al daño progresivo de las fibrillas que componen las fibras de colágeno. El puente entre la escala micro y macro, se establece a través de una técnica de homogeneización computacional con una integración numérica a lo largo de la superficie de la esfera de radio unidad. Dicho modelo se ha implementado en un código de elementos finitos a través de una subrutina de usuario de Abaqus. Se han ajustado los parámetros del modelo con ensayos experimentales y se han simulado diferentes casos y geometrías para comprobar el comportamiento del modelo. Por último, se ha simulado un proceso de angioplastia en una geometría real, donde los parámetros del modelo han sido ajustados con ensayos experimentales. 1 Contents 1 Introduction 5 2 Basic theory framework 7 2.1 Essentialkinematics .......................... 7 2.2 Hyperelastic framework . . . . . . . . . . . . . . . . . . . . . . . . . 8 3 Material model 11 3.1 Micro-sphere based model . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 Inclusion of anisotropy . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.2.1 The von Mises distribution . . . . . . . . . . . . . . . . . . . 13 3.2.2 The Bingham distribution . . . . . . . . . . . . . . . . . . . 13 3.3 Materialbehavior............................ 14 3.3.1 Mechanical behavior of the micro-fibers . . . . . . . . . . . . 15 3.4 Fitting material parameters and comparison . . . . . . . . . . . . . 16 3.5 Damagemodel ............................. 23 4 Validation examples 27 4.1 A displacement driven brick example . . . . . . . . . . . . . . . . . 27 4.2 Thin perforated plate . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5 Clinical application example: Simulation of angioplasty 33 6 Conclusions 37 A Obtention of stress and elasticity tensors 41 Bibliography 45 3 1 Introduction Most biological soft tissues, and particularly blood vessels, are composed of networks of collagen fibrils bundles (Rhodin, 1980) embedded in an isotropic ground substance with a high water content, which provide a well known almost incompressible behavior (Carew et al., 1968; Chuong & Fung, 1984). More in dept, blood vessels have three main layers (see, e.g., (Fung, 1990)), intima, media and adventitia. Media is mainly composed of smooth muscle and sheets of collagen fibers oriented preferentially along the circumferential direction. The adventitia, composed basically by a more random distribution of collagen fibrils bundles and the intima is made up of a thin layer of endothelial cells. The preferred orientation of the collagen fibers are main responsible of the anisotropic response and a highly non-linear behavior as detailed in numerous works (see, e.g., (Fung, 1990), (Humphrey, 1995)). It is worth noting the random distribution within the differentiable orientation of the collagen fibers in each layer as mentioned above. This structure on the micro-level has an important relevance in order to characterize these materials as will be discussed later. Many constitutive models have been proposed over the last years to characterize biological tissues, to model them in incompressible or quasi-incompressible hyperelastic frameworks. In this context, the description of a given material comes from the definition of a strain energy function (SEF) from which all the mechanical relations and variables can be obtained. Early SEFs in soft tissue mechanics were purely phenomenological, whereby non micro-structural information is gathered, so they just were able to describe some aspects of the material and, usually, only under physiological loads. However, many of them have been widely used to this purpose with really satisfactory results (Demiray et al., 1988). In order to provide a more realistic characterization of the tissue, structural models, in which a structural tensor is introduced in the free energy function to take into account the anisotropy of the material (see, e.g., (Boehler, 1987; Holzapfel et al., 2005; Menzel & Steinmann, 2003)), were proposed to deal with these limitations. Although it was an advance on the mechanical characterization of materials, and soft tissue in particular, many of the latest works about these topics are going toward a more micro structural characterization and therefore micro-structural 5 6Introduction models. Probably due to the improvements on the experimental field which have led to get information to feed these models. For example, the works of Ateshian (2007) and Taber (1998) are related to growth of the micro constituents, Kuhl & Holzapfel (2007) and Himpel et al. (2008) presented some results on the remodeling of collagen and elastin. Moreover, Gasser et al. (2006) and later Menzel et al. (2007) included micro-structure information into the hyperelastic formulation through the assumption of a statistical distribution of the fiber orientation around a preferential direction. In short, the high complexity of biological tissues requires mechanical models that include information of the underlying constituents and look for the physics of the whole processes within the material. The behavior of the micro-constituents can be taken into the macroscopic models by means of computational homogenization. It is in this context where the micro-sphere-based approach takes on a higher relevance. Miehe et al. (2004) used the micro-sphere approach with emphasis on elastomer and Alastrué et al. (2009a) carried it out applied to biological tissue. Usually, the use of these models to characterize this kind of tissues are limited to the range of physiological loads or elastic regions. Most of the finite element implementations carried out up today have been limited to these regions and just few of them deal with the failure of soft tissues. Damage over biological tissues is an important issue when dealing with these materials. Some physiological processes, as aneurysm or over-stretching of tendons, or clinical surgeries as angioplasties or clamping could lead to damage of the material. Damage in soft biological tissue could be produced by the progressive failure or softening of either the matrix or the fibers. In this way some authors have proposed damage models for only one component of the material (fibers or matrix), applied to transversal anisotropy of soft tissue (Natali et al., 2005; Ehret & Itskov, 2009; Peña & Doblaré, 2009). In this work, a three-dimensional decoupled finite strain formulation with a multi-scale model for the anisotropic part incorporating continuum damage is presented. The damage model is introduced both in the matrix and the collagen fibers. Though the damage along the fibers have been showed the most main issue, the damage over the matrix can not be avoided for a real modeling of the material. With the micro-structural incorporation a more realistic response of the damage evolution is hoped due to a transition of the damage in the micro scale from fibril to fibril all over the micro-sphere surface. The work is organized as follows. In Chapter 2 the basis of the mechanical framework is presented. Chapter 3 deal with the material model, including the micro-sphere approach, the mechanical behavior, the anisotropy and the damage models. Chapter 4 present some numerical examples with some sensitivity analysis of parameters in different geometries. Chapter 5 shows the results of a clinical interest example, and finally, Chapter 6 present the conclusions of this work. 2 Basic theory framework 2.1 Essential kinematics Let Ω0be the reference or material configuration of a continuum body Band Ω the current or spatial configuration at time t, regarding an arbitrary reference system (Cartesian since now on). Let X∈Ω0be the position of a particle in the reference configuration and x∈Ωthe position of the same particle at time t. The non-linear application ϕ:X→xthat relates the position of xto X,x=ϕ(X) is called motion map. The gradient of ϕrespect to X,F=∇Xϕ, represents the deformation gradient bi-point tensor. The Jacobian of the motion, J, is given by J=det(F). In order to reproduce the quasi-incompressible behavior of soft tissues, an uncoupled representation of the SEF is used (Flory, 1961). The multiplicative decomposition of the deformation gradient and the Cauchy-Green tensor C=FT·F can be expressed as F= [J1/3I]·F(2.1) C= [J2/3I]·C(2.2) where the terms J1/3Iand J2/3I, with Ithe second order unit tensor, are associated with changes of volume while Fand C, the isochoric deformation gradient and Cauchy-Green tensors, respectively, representing the volume-preserving part. Furthermore, let r∈Ω0be a vector of the reference configuration. The motion map turns this vector, with the so called push-forward operation, into overlinet∈ Ω. Assuming ris affected just by the isochoric part of F ¯ t=F·r=J−1/3twith k¯ tk=λ=J−1/3ktk(2.3) ¯ tcorresponds to the isochoric push-forward of the material vector rand ¯ λto the isochoric stretch in the direction of r. 7 14 Material model where Ais a symmetric 3x3matrix, dAis the Lebesgue invariant measure on the unit sphere, r∈U2and K(A)is a normalizing constant. As its main features, it is worth noting that this distribution always exhibits antipodal symmetry, but not rotational symmetry for the general case (Bingham, 1974). Applying straightforward transformations, Eq. (3.8) can be rewritten as ρ(r;Z,Q)dA 4π= [F000(Z)]−1etr A·At·r·rt·QdA 4π,(3.9) where etr (•)≡exp (tr (•)),Zis a diagonal matrix with eigenvalues κ1,2,3,Q∈Q3 such that A=Q·Z·QTand F000(Z)may be written as F000(Z) = [4 π]−1ˆU2 etr Z·r·rtdA=1F1(1 2;3 2;Z),(3.10) with 1F1a confluent hyper-geometric function of matrix argument as defined by Herz (1955). Thus, the probability concentration is controlled by the eigenvalues of Z, which might be interpreted as concentration parameters. Specifically, the difference between pairs of κ1,2,3– i.e., [κ1−κ2],[κ1−κ3]and [κ2−κ3]– determines the shape of the distribution over the surface of the unit sphere. Therefore, the value of one of these three parameters may be fixed to a constant value without reducing the versatility of (3.9). Figure 3.2 shows different distributions of the fiber bundles concentration achieved for a constant value of κ1and varying values of κ2and κ3. With this assumption, the distribution can expand along the plane composed of the associated direction with κ2and κ3. As seen in Fig. 3.2(a) and (b), setting two of the parameters equal to zero the Von Mises ODF is obtained. The highest value of the parameters lead to the preferential orientations as shows Fig. 3.2(c) and (f). As two parameters come close up to a value a rotational symmetry is achieved. 3.3 Material behavior As discussed above, the definition of a given material in the hyperelasctic modeling framework lies in the setting of free energy functions associated to each part of the above discussed splitting. Ψ = Ψvol(J)+Ψiso(I1)+Ψani(n, ρ, ¯ λ)(3.11) Ψvol(J) = 1 Dln2(J)(3.12) Ψiso(I1) = µ[I1−3] (3.13) Ψ(n, ρ, ¯ λ) = hnρfψfi(3.14) The ground substance is known to be composed by a high water content, which results in an almost incompressible behavior, so the volumetric part of the energy function supposed force the quasi-incompressibility depending on the value Material behavior 15 (a) κ2= 5.0,κ3= 5.0(b) κ2= 0.0,κ3= 50.0(c) κ2= 40.0,κ3= 50.0 (d) κ2= 49.0,κ3= 50.0 (e) κ2= 10.0,κ3= 10.0(f) κ2= 50.0,κ3= 49.0 Figure 3.2: Representation of the Bingham ODF for different sets of parameters with κ1= 0.0. of the penalty parameter D (3.12). The matrix contributes to the overall behavior through the volumetric and the isotropic part of the energy function (3.13). Regarding the anisotropic part of the model (3.14), two statistical dispersion functions of the fibers, around a preferential orientation, are considered . 3.3.1 Mechanical behavior of the micro-fibers Two microscopic strain energy density functions were used to model the contribution of the collagen fibrils to the macroscopic mechanical quantities (Alastrué et al., 2009a). First, a phenomenological exponential function frequently used to model the fiber response from a macroscopic phenomenological approach (Holzapfel et al., 2000), namely nΨi f(¯ λi) =    0if ¯ λi<1 k1 2k2hexp k2¯ λ2 i−12−1iif ¯ λi≥1(3.15) where k1is a stress-valued constant, k2is a dimensionless parameter, and λidenotes the isochoric stretch in the fiber direction of ri, i.e., λi=ktik. The particularization of the eight-chain model (Arruda & Boyce, 1993) to the transversely isotropic case (Kuhl et al., 2006; Alastrué et al., 2009a) was also used 16 Material model to model microfiber contributions, i.e.: nΨi f(¯ λi) =                  0if ¯ λi<1 n KΘ L 4A2¯r2 i L2+1 1−¯ri/L −¯ri L −ln(¯ λ4r2 0 i) 4r0Lh4r0 L+1 [1 −r0/L]2−1i−Ψrif ¯ λi≥1 (3.16) with ¯ri=¯ λir0, and Ψr= 2 r2 0 L2+1 1−r0/L −r0 L(3.17) being a repository constant accounting for a zero strain energy at ¯ λi= 1 (Alastrué et al., 2007). Notice that in both cases fibers are assumed not to bear any load under compression. 3.4 Fitting material parameters and comparison The strain energy functions in Section (3.3) was used to fit experimental curves of simple tension tests carried out on human coronary arteries (Holzapfel et al., 2005). Data coming from thirteen individuals for media and adventitia layers were fit. To be specific, experiments on samples along the circumferential and longitudinal directions of the cylindrical frame {eθ,z,r}were considered. A least-square scheme was used to minimise the functional χ2= p X j=1 hσθθ −σΨ θθ2 j+σzz −σΨ zz2 ji,(3.18) with pthe number of experimental measures, σθθ and σzz the Cauchy stress data obtained from the tests, and σΨ θθ and σΨ zz representing the Cauchy stresses obtained from (3.35) as σΨ=J−1τΨ. Nevertheless, in the case of the Bingham distribution, a single ODF was used to reproduce the curves where two are needed in the von Mises distribution, for which the eigenvalues of the Zmatrix – namely, κ1,2,3– were assumed to be in the interval [0,∞)representing the concentration parameters associated to the radial, the longitudinal and the circumferential directions, respectively. Then, following experimental evidences (Rhodin, 1980), the preferential orientation of the microfiber were assumed to be contained in the plane formed by vectors eθ and ezby fixing κ1= 0. A quasi-Newton minimization algorithm was used to minimize Eq. (3.18). The normalized root mean square error ε=1 µpχ2/[p−q],(3.19) Fitting material parameters and comparison 17 1 1.1 1.2 1.3 1.4 0 50 100 150 200 λ σ (KPa) Averaged constants Experimental data (a) Adventitia: longitudinal direction. 1 1.1 1.2 1.3 1.4 0 50 100 150 200 λ σ (KPa) Averaged constants Experimental data (b) Adventitia: circumferential direction. 1 1.1 1.2 1.3 1.4 0 50 100 150 200 λ σ (KPa) Averaged constants Experimental data (c) Media: longitudinal direction. 1 1.1 1.2 1.3 1.4 0 50 100 150 200 λ σ (KPa) Averaged constants Experimental data (d) Media: circumferential direction. Figure 3.3: Exponential fiber behavior: simulation results of uniaxial tension tests and experimental data. was used as a measure of the quality of the approximated set of parameters, where qis the number of parameters to be identified so that p−qis the number of degrees of freedom, and µrepresents the mean stress defined as µ=1 p p X j=1 [σθθ +σzz]j.(3.20) Following Alastrué et al. (2009a), the deformation gradient tensors causing the stress-strain curves reported by Holzapfel et al. (2005) were assumed to correspond to pure shear deformations, namely F=λeθ⊗eθ+λ−1ez⊗ez+er⊗erand F=λ−1eθ⊗eθ+λez⊗ez+er⊗erfor the circumferential and longitudinal samples, respectively. The obtained sets of constant for each of the individuals together with the averaged values and the standard deviation are collected, for the von Mises distribution, in Table 3.1 and 3.2 for the exponential micro-structural behavior (3.15). In Table 3.3 and 3.4 for the worm-like structural behavior in Eq. (3.16). The graphic representation is presented in Fig. 3.3 and Fig. 3.4. The mean and the standard deviation for the Bingham distribution are presented, in the same way, in Table 3.6 and Table 3.5. Three-dimensional plots of the von Mises and Bingham ODF, for the identified parameters presented above are depicted in Figure 3.5 and Figure 3.6 respectively. Note the similar values and shapes of the ODFs rendered by the identified set of constants for both micro-structural behaviors. As it is seen in Fig. 3.5a and b and Fig. 3.6a and b, the fibers are preferentially oriented along the longitudinal direction in the adventitia layer, whereas its preferred orientation is along the circumferential direction in the medial layer, Fig. 3.5c and d and Fig. 3.6c and d. 18 Material model Table 3.1: Material parameters identified for the adventitia layer based on exponential fiber behavior for the von Mises distribution: µ[kPa], k1[kPa], k2[-], b[-], ϑ[◦], ε[-]. Specimen µ k1k2b ϑ ε 1 4.04 82.768 36.954 10.340 65.43 0.187 2 2.41 79.657 13.721 10.616 60.84 0.262 3 4.91 89.441 14.154 8.271 55.17 0.117 4 16.41 18.606 11.260 4.502 68.13 0.149 5 9.65 32.094 16.797 2.936 58.74 0.139 6 7.92 78.434 22.478 14.038 63.27 0.250 7 7.07 90.307 12.062 5.187 67.05 0.242 8 2.25 48.856 12.156 5.954 59.22 0.075 9 6.16 41.872 24.257 3.149 67.77 0.088 10 14.93 91.374 55.450 15.000 51.99 0.168 11 12.71 76.972 29.472 10.315 62.37 0.155 12 3.13 13.411 8.779 1.582 60.83 0.046 13 6.69 21.935 5.321 1.998 68.58 0.189 Mean 7.560 58.902 20.220 7.222 62.260 0.159 SD 4.658 30.071 13.795 4.528 5.165 0.067 Table 3.2: Material parameters identified for the medial layer based on exponential fiber behavior for the von Mises distribution: µ[kPa], k1[kPa], k2[-], b[-], ϑ[◦], ε[-]. Specimen µ k1k2b ϑ ε 1 0.94 10.541 4.486 1.589 23.10 0.042 2 1.09 16.071 3.309 1.587 28.38 0.083 3 1.79 12.863 3.917 1.059 25.30 0.118 4 1.63 21.501 4.214 2.588 28.16 0.047 5 2.54 15.245 4.506 1.157 18.70 0.075 6 0.73 23.962 3.203 1.316 16.83 0.049 7 0.93 17.825 3.598 0.768 11.07 0.061 8 0.32 26.119 2.983 0.777 11.44 0.031 9 2.31 5.301 5.408 2.564 27.39 0.031 10 0.85 25.223 1.910 0.435 18.04 0.059 11 1.21 7.210 3.754 1.772 26.73 0.087 12 0.96 20.640 3.218 0.642 22.55 0.032 13 1.19 19.013 4.602 2.596 24.42 0.059 Mean 1.268 17.040 3.778 1.451 21.700 0.056 SD 0.634 6.664 0.894 0.755 5.996 0.026 Fitting material parameters and comparison 19 Table 3.3: Material parameters identified for the adventitia layer based on wormlike-chain fiber behavior for the von Mises distribution: µ[kPa], B[kPa], r0[mm], L[mm], b[-], ϑ[◦], ε[-]. Specimen µ B r0L b ϑ ε 1 4.04 3.441 0.527 0.677 4.311 65.43 0.337 2 2.41 3.741 0.100 0.139 2.812 74.36 0.351 3 4.91 2.344 0.135 0.179 2.244 67.43 0.190 4 16.41 1.360 0.597 0.813 2.289 74.89 0.240 5 9.65 0.843 0.390 0.491 2.075 58.74 0.245 6 7.92 5.277 1.089 1.512 3.049 77.33 0.415 7 7.07 4.221 0.100 0.135 2.209 81.95 0.317 8 2.25 1.886 0.512 0.681 2.404 72.38 0.238 9 6.16 0.664 0.233 0.287 1.564 82.83 0.187 10 14.93 1.520 0.251 0.312 4.010 58.74 0.333 11 12.71 2.696 0.100 0.130 2.840 76.23 0.302 12 3.13 0.757 0.766 1.027 0.623 55.14 0.159 13 6.69 1.589 0.581 0.817 1.069 82.12 0.325 Mean 7.560 2.334 0.414 0.554 2.423 71.352 0.281 SD 4.658 1.455 0.302 0.419 1.033 9.458 0.076 Table 3.4: Material parameters identified for the medial layer based on wormlike-chain fiber behavior for the von Mises distribution: µ[kPa], B[kPa], r0[mm], L[mm], b[-], ϑ[◦], ε[-]. Specimen µ B r0L b ϑ ε 1 0.940 0.561 0.968 1.340 1.601 21.0 0.059 2 1.090 1.042 1.077 1.538 1.484 25.8 0.161 3 1.790 1.013 1.190 1.695 0.945 23.0 0.212 4 1.630 0.962 0.952 1.314 2.440 25.6 0.109 5 2.540 0.852 1.104 1.534 1.101 17.0 0.162 6 0.730 1.096 0.912 1.277 1.494 18.7 0.061 7 0.930 0.865 1.077 1.506 0.834 12.3 0.050 8 0.320 1.411 0.972 1.380 0.805 10.4 0.062 9 2.310 0.570 1.023 1.439 2.748 24.9 0.224 10 0.850 1.824 1.145 1.700 0.440 16.4 0.099 11 1.210 0.474 1.158 1.637 2.486 29.7 0.145 12 0.960 1.263 1.012 1.442 0.660 20.5 0.080 13 1.190 1.310 0.992 1.395 2.393 22.2 0.199 Mean 1.268 1.019 1.045 1.477 1.495 20.576 0.125 SD 0.634 0.378 0.087 0.139 0.787 5.544 0.062 Table 3.5: Material parameters identified based on exponential fiber behavior for the Bingham distribution. Layer - µ[kPa] k1[kPa] k2[-] κ2[-] κ3[-] ε[-] Adventitia Mean 7.560 55.176 14.755 63.534 57.334 0.186 SD 4.658 29.106 8.638 16.301 14.355 0.072 Media Mean 1.268 19.622 3.437 61.606 63.416 0.072 SD 0.634 7.867 0.810 19.420 19.635 0.029 20 Material model Table 3.6: Material parameters identified based on worm-like-chain fiber behavior for the Bingham distribution. Layer - µ[kPa] B[kPa] r0[mm] L[mm] κ2[-] κ3[-] ε[-] Adventitia Mean 7.560 3.622 0.752 1.000 61.210 57.830 0.273 - SD 4.658 2.529 0.425 0.564 15.142 15.505 0.079 Media Mean 1.268 1.406 1.072 1.532 56.809 58.699 0.117 - SD 0.634 0.527 0.496 0.705 17.311 17.911 0.058 1 1.1 1.2 1.3 1.4 0 50 100 150 200 λ σ (KPa) Averaged constants Experimental data (a) Adventitia: longitudinal direction. 1 1.1 1.2 1.3 1.4 0 50 100 150 200 λ σ (KPa) Averaged constants Experimental data (b) Adventitia: circumferential direction. 1 1.1 1.2 1.3 1.4 0 50 100 150 200 λ σ (KPa) Averaged constants Experimental data (c) Media: longitudinal direction. 1 1.1 1.2 1.3 1.4 0 50 100 150 200 λ σ (KPa) Averaged constants Experimental data (d) Media: circumferential direction. Figure 3.4: Worm-like chain fiber behavior: simulation results of uniaxial tension tests and experimental data. Fitting material parameters and comparison 21 (a) Total ODF for the adventitia, exponential model: ϑ= 62.26◦, b = 7.606. (b) Total ODF for the adventitia wormlike-chain model: ϑ= 71.352◦, b = 2.423. (c) Total ODF for the media, exponential model: ϑ= 21.700◦, b = 1.451. (d) Total ODF for the media, worm-like- chain model: ϑ= 20.576◦, b = 1.495. Figure 3.5: Representation of the von Mises ODF, ρtot =ρ1+ρ2, for the identified averaged values band ϑ– the angle ϑfor the purpose of visualization being referred to the vertical direction. 6 eθ H Hj ez   er The von Mises and the Bingham statistical functions has been considered in this work for the incorporation of anisotropy to a hyperelastic micro-sphere-based constitutive law with application to the modeling of the vascular tissue. The obtained results show a good agreement with the experimental data, with values of the normalized root means squared error very similar between them. Homogeneous deformation states were reproduced in order to compare both the von Mises and Bingham ODFs. An equivalent response is attached for the two ODFs for the uniaxial tests (Fig. 3.7(a) and 3.7(b) and Fig. 3.8(a) and 3.8(b)). For the equi-biaxial tests, a good correlations also found in all cases except for the exponential fiber behavior in the adventitia layer Fig. ( 3.7(c) and (d) and Fig. 3.8(c) and (d)). In view of these results, it is worth noting that the use of the Bingham distribution allowed getting similar results with a single family of fibers compared with the two helically oriented families of fibers – or rather preferential orientations – modeled by means of the von Mises ODF. This fact, which is crucial in the reduction of the computation associated the numerical calculation of the macroscopic stress, points out the feasibility of the Bingham distribution to include anisotropy in the micro-sphere based models. Nevertheless, some simplifications have been taken for the obtaining of the results here presented. The most important one regards the use of the eigenvalues 22 Material model (a) Adventitia, exponential model: κ2= 63.53,κ3= 57.33. (b) Adventitia, wormlike-chain model: κ2= 61.21,κ3= 57.83. (c) Media, exponential model: κ2= 61.61, κ3= 63.42. (d) Media, worm-like- chain model: κ2= 56.81,κ3= 58.70. Figure 3.6: Representation of the Bingham ODF for the identified averaged values κ2and κ3– the circumferential direction is referred to the ezdirection, whereas excorresponds to the longitudinal one. 6 eθ H Hj ez   er (a) Uniaxial test, media layer. (b) Uniaxial test, adventitia layer. (c) Biaxial test, media layer. (d) Biaxial test, adventitia layer. Figure 3.7: Comparison between Von Misses and Bingham ODF for exponential fiber behaviour (Media and adventitia layer). Damage model 23 (a) Uniaxial test, media layer (b) Uniaxial test, adventitia layer (c) Biaxial test, media layer (d) Biaxial test, adventitia layer Figure 3.8: Comparison between Von Misses and Bingham ODF for eight-chain model fiber behaviour (media and adventitia layer). of Zas additional identification parameters for macroscopic stress-stretch curves. On the contrary, those parameters should by identified from measurements at the micro-structural level (Landuyt, 2006). Moreover, the eigenvectors of Z were assumed coincident to the vectors eθ,eθand ez, when the actual orientation at the micro-structural level can follow a more random distribution. However, this simplification has a strong experimental basis, since most of the histological studies show that most of the fibers are comprised in that plane (Rhodin, 1980). Moreover, the contribution of the isotropic neo-Hookean term constant µwas fixed to value obtained by Holzapfel et al. (2005), which might influence the results of the fittings. Nevertheless, it has been reported that the isotropic contribution to the macroscopic stress is only relevant up to certain value of stretch, from which collagen fibers bear most of the load (Holzapfel et al., 2000). In addition, fixing µallowed us to establish direct comparisons with the values reported in Alastrué et al. (2009a). Finally, the deformation gradient causing the stress-stretch curves was assumed to correspond to a uniaxial and equibiaxial tension with incompressibility constraint. Other deformation fields fulfilling this requirement might have been used, though previous studies show that the obtained constants are barely dependent on the used deformation field (Alastrué et al., 2009a). 3.5 Damage model The structural model used in this work is able to describe the softening of the material in a large strain non-lineal framework using the concept of internal variables 30 Validation examples (a) Map of stress in X axis at 36.9%. (b) g at 50% (c) Sxat 36.9% (d) Map of stress in X axis at 63.8%. (e) [g at 63.8% (f) Sxat 63.8% Figure 4.3: Evolution of stress and [1 − D]in the integration point for case 1. (a) Map of stress in X axis at 50%. (b) g at 50% (c) Sxat 50% (d) Map of stress in X axis at 100%. (e) g at 100% (f) Sxat 100% Figure 4.4: Evolution of stress and [1 − D]in the integration point for case 2. Thin perforated plate 31 Results present some important differences in both cases. Case 1 presents a tougher response, since the damage starts later than Case 2, and can not achieve total damage. At 50% of strain Case 1 present a maximum integrity factor of 0.84 and at the end reaches 0.6 in just the integration directions more alignment with exaxis (note that when gincrease, the damage D decrease). Case 2, with a sooner starting damage showed a value of 0.1 at the 50% and with a spreader distribution and damage close to 1 is obtained at 100% of strain. It is worth noting the evolution of the stress over the stereographic projection. Fig. 4.4.b capture a situation where the fibers along the exaxis present a higher value than the other ones (although lower than those achieved without damage). However, Fig. 4.4.e shows a circular crown-like shape caused by the total failure of those fibrils oriented along exaxis. Moreover, the more strain is achieved the more movement of the crown towards the middle plane. 5 Clinical application example: Simulation of angioplasty One of the aims of these models is the realistic simulation of clinical applications. Angioplasty is one of the most widely spread techniques in vascular surgery. Therefore, the study of the mechanical behavior of vessels in this situations is an important task, since this procedure can affect to the vessel integrity reducing its stiffness (Oktay, 1994). In order to assess this behavior, an angioplasty procedure was simulated using the above presented model. A pig aorta artery was experimentally tested to fit the elastic and damage material parameters which were used in the finite element simulations. In the following section the material parameters of the constitutive behavior and the damage model are fitted using experimental data (Peña et al., 2010). The vessel sample were cut comprising a longitudinal and a circumferential strips in order to perform uniaxial tests, and a least-square iterative method, as discussed in Section 3.4, is used for the identification process. The fitting procedure was carried out in two steps, the first one to fit the elastic part (before damage) and the second one for the damage part, when the softening phenomenon has started. The analytical stress expression above mentioned needs a deformation imposed to be calculated. In this approach have been considered incompressibility behavior with a deformation gradient tensor given by F=λeθ⊗eθ+λ−1ez⊗ez+er⊗er, considering strain in the direction of pulling (eθ),for the circumferential sample, or ezin longitudinal direction. Following these assumptions, Fig. 5.1 shows the identified parameter for the experimental test performed by Peña et al. (2010). The goal of this study is not the simulation of a rigorous clinical study with accurate. Instead of that, the objective is to show the applicability of the preµ[KPa] k1[KPa] k2[-] b[-] a[-] c[KPa] 1.051 246.99 0.9111.49 52.77 0.61 26.10 Table 5.1: Elastic and damage parameters for the uniaxial test. 33 34 Clinical application example: Simulation of angioplasty (a) Undeformed shape. (b) Deformed shape. Figure 5.1: Undeformed and deformed shapes for the angioplasty simulation. sented model to simulated the vessel behavior. The finite element simulation of the angioplasty was performed in a eight part of the model, applying symmetry conditions. The model consist of a 10 mm length sample with an external diameter De= 5mm and an internal diameter Di= 3.7mm corresponding to coronary artery dimension. One only layer was considered since the well know intima, media and adventitia layers were not divided in the experimental test, so only information of the damage behavior of the whole artery was gathered. The artery is composed by an isotropic ground substance and two families of fibers, the angle respect to the circumferential direction show in Fig. 5.1 and the parameters shown in the same table. The geometry is discretized in 4788 hexahedrals elements. The balloon present the same dimension and material properties than in Alastrué et al. (2007); Gasser & Holzapfel (2007). The step loads are applied sequentially as follows: (i) Imposition of a initial deformation gradient as proposed by Rodriguez et al. (1994); Alastrué et al. (2007), (ii) Application of internal pressure of 13.3 kPa in the vessel supposing as the mean physiological pressure condition and (iii) Imposition of pressure to the internal face of the balloon following the curve shown in Alastrué et al. (2007) in order to reach the contact between vessel and balloon. The undeformed and deformed shapes of the model are presented in Fig. 5.1a and 5.1b respectively. The maximal principal stress map on the artery at the end of the analysis is presented in Fig. 5.2a. The maximum stress, due to the one only layer simplification, is placed on the inner radius of the vessel. Some other authors (Alastrué et al., 2007), due to the incorporation of two different layers, observed different distribution of stress. Since the adventitia layer have a more circumferential distribution of the fiber, the maximal principal stress had higher values than the media layer. In order to present a mean value of the damage distribution over a given integration point depicted in Fig. 5.2b an averaged or 35 (a) Maximal principal stress map. (b) Integrated damage map. Figure 5.2: State of the vessel at the end of the analysis integrated value of the damage, defined as Dave =ˆU2 ρDdA ≈ m X i=1 ρiDiwi(5.1) According to the stress results, the highest damage value is localized on the inner radius. For a better visualization of the behavior of the vessel, Fig. 5.3 presents the the micro-structure stress at the end of the analysis for the undamaged case, Fig. 5.3a, and the damaged case, Fig. 5.3b. It is clearly appreciable the two family of fibers stress and the influence of the damage phenomenon. Fig. 5.4 presents the stress in the circumferential direction, showing that the vessel start to damage in the last increments of the analysis , observing a softening effect. (a) Micro stress for the undamage vessel (b) Micro stresses for the damaged vessel Figure 5.3: Micro-stresses for the undamaged and damaged model. 36 Clinical application example: Simulation of angioplasty Figure 5.4: Maximum stress for circunferential direction for the damaged and undamaged model. 6 Conclusions The main goal of this work was to present a formulation for a damage model within a micro-sphere-based approach in order to a better characterization of this phenomenon in biological soft tissues. The model was formulated in a hyperelastic framework, using the splitting of the free Helmholtz energy density function into a volumetric and an isochoric part, which was split again into a isotropic part, associated with the extracellular matrix and an anisotropic part which meanly deals with the distribution of the collagen fibers. To this end, the free energy associated to the fibrils was expressed in terms of the stretch rate instead of the widely used invariant form. Therefore, in the framework of a multi-scale model an homogenization scheme was required to move from the micro to the macro level. The micro-sphere-based approach was used to this purpose, carrying out a numerical integration of the stress and elastic tensors over the surface of the unit sphere. Moreover, within this approach, a phenomenological damage model, previously developed for unidirectional families of fibers was incorporated in the micro-sphere-based approach. Previous works (Alastrué et al., 2009a) studied different discretization for the numerical integration, leading to 368 discretization as a optimum value in order to avoid a lack of accuracy in the results. The model was implemented in a commercial finite element code (ABAQUS) through a material subroutine (UMAT) and checked with different sets of material parameters and geometries, in particular a brick and a perforated plate, since it is well known that damage models present convergence problems, in particular in driven force problems. In fact, for the uniform strain example, convergence was achieved for all the sets of parameters, while it was not the case for the plate due to localization problems. The results reported in the brick case showed some important characteristics of the proposed model. As far as the damage parameters concern, the parameter a(related on how fast the evolution of the damage is) showed that the lower the parameter the smoother the evolution on the stress; while c(energy level at which the damage starts) showed that the lower the value of cthe lower the stress achieved. The higher the concentration parameter the higher the maximum value of stress and the lower the point of total failure. This point will be discuss later. In order to obtain the behavior in a non homogeneous deformation 37 38 Conclusions state in the thin perforated plate, other two sets of parameters was investigated. Classical convergence problems due to localization and loss of ellipticity appear and therefore, depending on the chosen parameters the convergence rate increases or decreases. Apart from this issue, the result showed high level of damage and stress concentration around the hole. Regarding the angioplasty simulation, in spite of lack in experimental data and the mono-layer simplification, it demonstrates the capability of the model to simulate clinical applications, including damage phenomenon, with higher level of information. Although the fitting scheme is still a handicap and some numerical integration improvements could be done in the future, these micro-structure models present a chance in order to know what is happening inside the material. Therefore, the present model has some important limitations. The most important one is the undetermined parameter fitting problem. The idea followed in this work to fit the experimental data with a least-square minimization algorithm being aware that the achieved parameter sets are not the unique ones able to reproduce the material behavior. Although some restrictions were applied in order to reduce the solution field, the authors couldn’t get a unique solution. The best way to reduce this underdetermination would be a direct measurement of the micro-structural parameters, through experimental tests. Other limitation, which avoid an ideal development of the model is the problem coming from the mismatch between the distribution of the concentration parameter and the integration directions. In cases with high concentration distributions, there are many of the integration direction that are not used in a optimum way, since match whit position where no fibrils exits. Alastrué et al. (2009b) worked about this issue and will be a priority line of future work. The last important problem arises from the well know lost of ellipticity of the problem on damage simulations. Developing of a model that overcomes these localization problems (Steinmann, 1999) is left for future work. In spite of these limitations, the one dimensional character of the constitutive equation applied at the micro-level offers huge possibilities, due to the simplicity and the possibility of incorporation of other micro-variables. Menzel & Waffenschmidt (2009) have reported some works related to remodeling processes and the works Göktepe & Miehe (2005); Miehe & Göktepe (2005) developed some inelastic models on the micro-sphere framework for isotropic materials. Moreover, the incorporation of the von Mises or the Bingham ODF could allow, in a really physical way, the development of some remodeling models such as those reported by Kuhl et al. (2005) or Menzel et al. (2007), whose preferential orientation direction can evolve during the simulation. This probabilistic function may be used coupled with growing and remodeling models, accounting for the mass transference and also reorientation of the fibers that could be modeled, for example, by modification of the ODF of the fibers. Therefore, it seems clear that the research field including these micro structure models is extremely wide and can be improve the constitutive models of soft biological tissue. As far as the comparison between the two ODF studied, the presented results show the capability of the presented micro-sphere models together with the Bing- 39 ham or the von Mises distribution to reproduce the mechanical behavior of soft biological tissues, in particular, of the human coronary artery. In conclusion, the model presents a multi scale model capable to reproduce the softening behavior of biological soft tissues in general and blood vessels in particular. The present model uses a micro-structural approach, which represent a clear advantage over other phenomenological damage models, since almost all of them treat fibers as a unidimensional material that can not capture the progressive failure of the fibrils bundle. C. J. Chuong & Y. C. Fung (1984). ‘Compressibility and constitutive equation of arterial wall in radial compression experiments’. J Biomech 17(1):35–40. H. Demiray, et al. (1988). ‘A stress-strain relation for a rat abdominal aorta’. J Biomech 21(5):369–374. A. E. Ehret & M. Itskov (2009). ‘Modeling of anisotropic softening phenomena: Application to soft biological tissues’. Int J Plasticity 25(5):901–919. P. J. Flory (1961). ‘Thermodynamic relations for high elastic materials’. T Faraday Soc 57:829–838. Y. C. Fung (1990). Biomechanics: Mechanical Properties of Living Tissues. Springer. T. C. Gasser & G. A. Holzapfel (2007). ‘Finite element modeling of balloon angioplasty by considering overstretch of remnant non-diseased tissues in lesions’. Comput Mech 40(1):47–60. T. C. Gasser, et al. (2006). ‘Hyperelastic modelling of arterial layers with distributed collagen fibre orientations’. J Roy Soc Interface 3:15–35. S. Göktepe & C. Miehe (2005). ‘A micro-macro approach to rubber-like materials. Part III: The micro-sphere model of anisotropic Mullins-type damage’. J Mech Phys Solids 53(10):2259–2283. C. S. Herz (1955). ‘Bessel Functions of Matrix Argument’. Ann Math 61(3):474– 523. G. Himpel, et al. (2008). ‘Time-dependent fibre reorientation of transversely isotropic continua - Finite element formulation and consistent linearization’. Intl J Numer Meth Eng 73(10):1413–1433. G. A. Holzapfel (2000). Nonlinear Solid Mechanics: A Continuum Approach for Engineering. John Wiley & Sons. G. A. Holzapfel, et al. (2000). ‘A New Constitutive Framework for Arterial Wall Mechanics and a Comparative Study of Material Models’. J Elasticity V61(1):1– 48. G. A. Holzapfel, et al. (2005). ‘Determination of layer-specific mechanical properties of human coronary arteries with nonatherosclerotic intimal thickening and related constitutive modeling’. Am J Physiol Heart Circ Physiol 289(5):H2048– 2058. J. D. Humphrey (1995). ‘Mechanics of the arterial wall: Review and directions’. Crit. Rev. Bio. Eng. 23(1-2):1–162. 46 E. Kuhl, et al. (2005). ‘Remodeling of biological tissue: Mechanically induced reorientation of a transversely isotropic chain network’. J Mech Phys Solids 53(7):1552–1573. E. Kuhl & G. Holzapfel (2007). ‘A continuum model for remodeling in living structures’. J Mater Sci 42(21):8811–8823. E. Kuhl, et al. (2006). ‘On the convexity of transversely isotropic chain network models’. Philos Mag 86:3241–3258. M. Landuyt (2006). ‘Structural quantification of collagen fibers in Abdominal Aortic Aneurysms’. Master’s thesis, Royal Institute of Technology in Stockholm, Department of Solid Mechanics and Ghent University, Department of Civil Engineering. A. Menzel, et al. (2007). ‘Towards an orientation–distribution–based multi–scale approach for remodelling biological tissues’. Computer Methods in Biomechanics and Biomedical Engineering . Accepted for publication. A. Menzel & P. Steinmann (2003). ‘A view on anisotropic finite hyper-elasticity’. Eur J Mech A/Solids 22(1):71–87. A. Menzel & T. Waffenschmidt (2009). ‘A microsphere-based remodelling formulation for anisotropic biological tissues’. Philos. Trans. R. Soc. London, Ser. A 367(1902):3499–3523. C. Miehe (1995). ‘Discontinuous and Continuous Damage Evolution In Ogden-type Large-strain Elastic-materials’. Eur. J. Mech. A. Solids 14(5):697–720. C. Miehe & S. Göktepe (2005). ‘A micro-macro approach to rubber-like materials. Part II: The micro-sphere model of finite rubber viscoelasticity’. J Mech Phys Solids 53(10):2231–2258. C. Miehe, et al. (2004). ‘A micro-macro approach to rubber-like materials–Part I: the non-affine micro-sphere model of rubber elasticity’. J Mech Phys Solids 52(11):2617–2660. A. N. Natali, et al. (2005). ‘Anisotropic elasto-damage constitutive model for the biomechanical analysis of tendons’. Med. Eng. Phys. 27(3):209–214. E. A. D. Neto, et al. (1998). ‘Continuum modelling and numerical simulation of material damage at finite strains’. Arch. Comput. Meth. Eng. 5(4):311–384. R. W. Ogden (1996). Non-Linear Elastic Deformations. Dover Publications. H. Oktay (1994). ‘Continuum damage mechanics of ballon angioplasty.’. In Internal Report. UMI. E. Peña, et al. (2010). ‘A constitutive formulation of vascular tissue mechanics including viscoelasticity and softening behaviour’. J Biomech 43(5):984–989. 47 E. Peña & M. Doblaré (2009). ‘An anisotropic pseudo-elastic approach for modelling Mullins effect in fibrous biological materials’. Mech Res Commun 36(7):784–790. J. A. G. Rhodin (1980). Handbook of Physiology, The Cardiovascular System, vol. 2, chap. Architecture of the vessel wall, pp. 1–31. American Physiological Society, Bethesda, Maryland. E. K. Rodriguez, et al. (1994). ‘Stress-dependent finite growth in soft elastic tissues’. J Biomech 27(4):455–467. J. C. Simo (1987). ‘On a fully three-dimensional finite-strain viscoelastic damage model: Formulation and computational aspects’. Comput Method Appl M 60(2):153–173. J. C. Simo & T. J. R. Hughes (1998). Computational Inelasticity. Springer. A. J. M. Spencer (1954). ‘Theory of Invariants’. In Continuum Physiscs, pp. 239–253. Academic Press, New York. P. Steinmann (1999). ‘Formulation and computation of geometrically non-linear gradient damage’. Intl J Numer Meth Eng 46(5):757–779. L. A. Taber (1998). ‘A model for aortic growth based on fluid shear and fiber stresses’. J Biomech Eng-T ASME 120(3):348–354. C. Truesdell & W. Noll (2004). The Non-Linear Field Theories of Mechanics. Springer-Verlag, 3rd edn. E. W. Weisstein (2004). ‘“Erfi.” From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/Erfi.html’ . 48