1
Optimizing numerical Aneurysm growth predictions
Gustav Thorén
Degree project in Solid Mechanics Second level, 30.0 HEC Stockholm, Sweden 2015
Abstract
Abdominal Aortic aneurysm (AAA) is a localized dilation of the abdominal aorta exceeding the normal diameter by more than 50 %. A common complication is rupture of the vessel wall, which in most cases is life-threatening. In order to prohibit loss of life, efforts in trying create computational models to predict failure of the vessel wall has been made. Current non-linear finite element models for biomechanical rupture risk assessment (BRRA) are very computational demanding and require very long computational times.
A BRRA using an automatic time-step scheme is presented in order to speed up the computational time of the very demanding Finite Element Analysis (FEA) models available. The numerical scheme is governed by the collagen dynamics and two numerical approaches are presented. One in which the time step depends only on the maximum incremental change of collagen in one of all directions, and a second in which the time step is governed by the total incremental change of collagen in all directions and the incremental change of the undulations stretches. Simulations were carried out both without and with contact conditions with the spine.
2
Optimering av
bukaortaaneurysmmodellering
Gustav Thorén
Examensarbete i Hållfasthetslära Avancerad nivå, 30 hp Stockholm, Sverige 2015Sammanfattning
Bukaortaaneurysm (AAA) är en lokal utvidgning av den abdominal aorta som är mer än 50 % större en den normala diametern. En vanlig komplikation är ruptur av artärväggen vilket oftast är livshotande. Försök till att förutspå risken för ruptur har gjorts, men de tillgängliga icke-linjära finita element metoderna för riskbedömning (biomechanical rupture risk assessment - BRRA) är mycket krävande och tar mycket lång tid att exekvera. En BRRA som använder en automatisk anpassning av tidssteget presenteras i denna uppsats som ett försök till att reducera beräkningstiden för en existerande BRRA. Tidsteget är beroende av den dynamiska kollagen-tillförseln. Två metoder presenteras, en där tidsteget bestäms av den maximala inkrementella ändringen av kollagen av alla riktningar. Den andra metoden summeras istället den inkrementella ändringen av kollagen i alla riktningar, dessutom tas kollagenets maximala och minimala sträckning i beaktning. Simuleringarna genomfördes både med och utan kontakt med ryggraden. Tidiga resultat visar en 50-75 % minskad beräkningstid med ett litet relativt fel. Simuleringarna resulterade också i ruptur vid spänningar som är fysiologiskt observerade. Detta tyder på att den konstitutiva modellen är sund och att tidstegsmetoderna har potential.
3
Examensarbete Solid Mechanics 2015-01-11
4
1. Nomenclature
Notations
Symbol
Description
E Young´s modulus (Pa)
t Time
Cauchy stress
T First Piola-Kirchoff stress
l Number of integration point
Stretch
Collagen density
CDF Cumulative density function
Elasticity tensor
LDE Local discretization error
, i er r o r
5
2. Table of Contents
1. Nomenclature ... 4 3. List of figures... 7 4. Introduction ... 8 4.1 Background ... 8 4.2 Purpose ... 8 4.3 Method ... 8 5. Frame of Reference ... 85.1 Introduction to an existing multiscale constitutive framework ... 8
5.1.1 Collagen model ... 9
5.1.2 Collagen dynamics ... 10
5.1.3 Macroscopic model ... 11
5.2 Numerical implementation ... 12
5.3 Automatic time stepping ... 14
5.3.1 Error estimation ... 14
5.3.2 Stability for explicit schemes ... 15
5.3.3 Examples of different automatic time stepping schemes ... 16
5.4 Basic theory Contact mechanics ... 17
5.4.1 Contact kinematics ... 17
5.4.2 Contact constraint treatments ... 18
5.5 Discretization of contact constraints ... 19
5.5.1 LAGRANGE multiplier ... 19
5.5.2 Penalty method ... 20
6. Method ... 20
6.1 Geometry and load case ... 20
6.1.1 One brick element with prescribed displacement ... 20
6.1.2 One brick element with prescribed load ... 21
6.1.3 Abdominal Aorta with internal pressure ... 21
6.1.4 Abdominal Aorta with internal pressure and contact with spine ... 21
6.2 Automatic time stepping ... 24
6.2.1 Automatic time stepping first approach ... 25
6.2.2 Automatic time stepping second approach ... 26
6.2.3 Estimation of error ... 26
6
7.1 First time stepping approach ... 27
7.1.1 Results: One brick element with prescribed displacement ... 28
7.1.2 Results: Abdominal Aorta with internal pressure ... 30
7.2 Second time stepping approach ... 32
7.2.1 Results: One brick element with prescribed load ... 32
7.2.2 Results: Abdominal aorta with internal pressure and contact with spine .... 33
8. Discussion and Conclusions ... 35
7
3. List of figures
Figure 1. Q1-P0 Element, for simplicity in 2D... 12
Figure 2. Contribution to average stresses in nodes. For simplicity in 2D. ... 13
Figure 3. Scheme for adaptive time stepping... 16
Figure 4. Contact illustration, Picture from (Kloosterman, 2002) ... 17
Figure 5. Parameterization of the contact surface, picture from (Kloosterman, 2002) 18 Figure 6. Geometry and load case 1: Q1-P0 brick ... 20
Figure 7. Geometry AAA ... 22
Figure 8: Proportional load, load case AAA with internal pressure. The load 13.6 kPa was constant until end of simulation. ... 22
Figure 9. CT scan of the AAA and spine ... 23
Figure 10. AAA and spine, von Mises stress. ... 23
Figure 11. Outline of the automatic time stepping scheme ... 25
Figure 12. Outline of the parameter test ... 26
Figure 13: Parameters for the first automatic time stepping scheme ... 28
Figure 14. Maximum principal stress and CPU-time, load case 1, constant time step .. 29
Figure 15. Relative error for different R-values, load case 1, compared to a constant time step of 0.001. ... 29
Figure 16. CPU-time for different R-values, load case 1. ... 30
Figure 17: Load case 2, results, element 802 ... 31
Figure 18. Relative error and CPU-time, AAA with internal pressure. ... 32
Figure 19. Time vs displacement, one brick element with prescribed load. ... 33
Figure 20. Time vs size of time-step, one brick element with prescribed load. ... 33
Figure 21. Size of time step vs time, AAA and spine, rhomax = 1.2. ... 34
Figure 22. Size of time-step vs time, AAA and spine, rhomax = 0.7. ... 34
Figure 23. Penetration of spine, showing von Mises stress [kPA]. ... 35
8
4. Introduction
4.1 Background
An Abdominal Aortic Aneurysm (AAA) is a localized dilation of the abdominal aorta exceeding the normal diameter by more than 50 % (Upchurch & Schaub, 2006). They occur most commonly in individuals between 65 and 75 years and are more common among men. Very few symptoms tend to be caused, they seldom cause pain in the abdomen and back. Most severe complications are ruptures, which are life-threatening - mortality of rupture is between 60 - 90 %.
In order to prohibit loss of life, efforts in trying create computational models to predict failure of the vessel wall has been made. The non-linear finite element models for biomechanical rupture risk assessment (BRRA) available are very computational demanding and takes very long time to run.
In order to speed up results from the BRRA, the time step should be adjusted to the solution process. As in all applications of computational models that involve rate equations, the error of the computed solution compared to the “real” solution depends on the size of the time step. Efforts are made in this thesis to connect the time step, and hence the error of the numerical computation, to a physical parameter. In the model described in this thesis it is related to the remodeling of the collagen in the AAA wall. Also, the impact of contact with the spine has not yet been incorporated in existing BRRA-models. The stress distribution in the AAA should be affected by the contact.
4.2 Purpose
The purpose is to implement an automatic time stepping algorithm and to test its validity. Also, the influence on the stress distribution in the AAA resulting from contact with the spine is investigated.
4.3 Method
Simulations are carried out using FEAP. Results are evaluated in Paraview and Microsoft Excel.
5. Frame of Reference
5.1 Introduction to an existing multiscale constitutive framework
The model builds on a modified multiscale constitutive framework presented by (Martufi & Gasser, 2012).
9 ,
vol i so
(5.1)
where i sois the isochoric part and volis the volumetric part. The volumetric part enforces the incompressibility, as known for vascular tissue, and is only depended of the volume dilation. The stored energy function of the vol is of a simple form
21/ 2 1 ,
vol J
(5.2)
where is the bulk modulus and J detF denotes the jacobian and F is the deformation gradient. For each time step is adjusted so that J 1 .
Specifically, the constitutive properties of the tissue are captured by iso, which in turn
defines its Cauchy stress . The tissue is regarded as a fibrous collagenous tissue; fibers of collagen reinforces an otherwise isotropic matrix material. The complete macroscopic Cauchy stress is composed of
,
vol nh f ibr e
(5.3)
where nh corresponds to the stress in the matrix material (Neo-Hookean) which is thought to capture the mechanical properties of mainly elastin and f ibr e corresponds to the stress in the collagen fibers.
5.1.1 Collagen model
The framework takes collagen turnover in consideration in order to accurately describe the macroscopically properties of AAAs.
Collagen fibrils are regarded virtually linear elastic and the undulation follows a triangular probability, which yields in turn an expression for the first Piola-Kirchoff stress according to
0 , T k CDF d
(5.4)where CDF
denotes the Cumulative Density Function of the triangular probability distribution and k denotes the stiffness of the proteoglycan cross-linked collagen fibrils (CFPG-complex), see (Gasser & Martufi, 2011).10
3 2 3 min 2 0 0 , 2 , 3 `, 2 , 3 , mi n mi n ma x ma x ma cfibr x e k k k (5.5)with
ma x
min and
ma xmi n
/ 2. The limits
min and
max are updated every iteration according to the active collagen dynamics.5.1.2 Collagen dynamics
The collagen fibre in the material is modeled by a collagen turnover model. As stress or strain rises fibroblast sense this and adjust the amount of collagen in the network. Collagen turnover is accomplished by fibroblast cells that are spread out throughout the collagen network and secreted by fibroblast and are assembled into organized suprafibrillar structures.
The fibroblasts sense the collagen fibre stretch . If one assumes that 10 % of the collagens fibrils are engaged in the AAA wall at physiological deformation and that the collagen fibrils are distributed in a triangular probability density function, the physiological stretch (the stretch at physiological deformation) is
/ 2 0.224 ,
min min
ph ph
(5.6)
where ph 0.1denotes the density of engaged collagen fibrils at the physiological deformation (Gasser & Martufi, 2011).
The condition of 10 % engaged collagen fibrils is regarded as the homoestatic condition. Consequently, newly formed collagen is deposited at
min 0.224 .
new new new
ph
(5.7)
In order to reach the stress new ph
collagen fibrils are deposited with
min 0.224
new new new
ph
and max 0.776
new new new
ph
where “new” denotes that these values are for the newly deposited collagen. Regarding a collagen fiber aligned with the referential direction N , an optimality condition arises at this homeostatis according to
11
where
N is a mechanical stimulus. For
N 1 , collagen is stretched too much, such that more collagen is required to reach homeostatis. Similarly for
N 1, a net loss of collagen is needed to reach homeostatis.The collagen production is assumed to be related to the mechanical stimulus according to
max min , , N (5.9) where max denotes the maximum collagen production rate- Here, is the total collagen density and denotes the time-scale of collagen production.
Collagen is assumed to degrade isotropically and independent of the stretch according to
,
(5.10)
such that the degradation is thus purely time-based. Finally, the net contribution is
,
(5.11)
and 0 zero at homeostasis.
The collagen fibrils are disintegrated from the current structure without changing the undulation characteristics. They do however, if the current slab is not at homeostasis, change the characteristics at integration. The fibrils where integrated at a certain distribution of pre-stretches. Further information is given in (Gasser & Martufi, 2011).
5.1.3 Macroscopic model
The macroscopic collagen stress in each Gauss point is computed by integrating over all fiber-directions (Gasser & Martufi, 2011). The integration of (5.12) is done by a spherical t-design. A t-design integrates a polynomial of exactly t degrees. Specifically the concept is that the spherical integration over the solid angle according to
,f ibr e dev d
N N nn (5.12)which can be approximated by
1 4 ... , int i n l l l t d l
(5.13)12
Table 1. Spherical t-design (J.Burkart, 2014)
Design Number of integration points in 3D
0 1
1 2
9 48
21 240
5.2 Numerical implementation
The framework is implemented using a P0 mixed finite element formulation. The Q1-P0 element is one of the mostly used element for incompressible or nearly incompressible materials.
The element is a hexahedral, with eight nodes and one single elemental pressure degree of freedom, see Figure 1. The element has eight displacement degree of freedoms and one pressure degree of freedom (Simo & Taylor, 1985).
Figure 1. Q1-P0 Element, for simplicity in 2D.
The macroscopic stress in the fibers, aggregated from all directions of the fibrils, is computed in every gauss point on the element. In 3d a 2x2x2 Gaussian integration order is implemented which results in eight Gauss points on each element, see Figure 2 The stress in each Gauss point is there after extrapolated using trilinear extrapolation to the nodes, taking in account the integration points in all the surrounding elements, see Figure 2.
13
Figure 2. Contribution to average stresses in nodes. For simplicity in 2D.
The outline of the numerical implementation of the calculation of the stresses at the Gauss point is presented below (Gasser & Martufi, 2011)
Initialization, t = 0
Spherical design is defined, i.e 0,1,9,21…-design, which governs the number of integration points lint, which governs the number of directions Ni where i1,...,li nt and the integration weight w4 /
lint.Also, the undulation stretches
min and ma x and the collagen fiber density is loaded. At time tRead the isochoric part of the deformation gradient F
Initialize the stress sensor of the fiber fiber and the elasticity tensor fibre Then a loop over all integration points is initialized:
DO i1,lint
Compute the fiber direction nFNi and n
Check if the undulation stretch is greater than the minimum undulation stress (i.e if the collagen fibril is in a load bearing mode)
IF min THEN
Compute the stress factor a
T / 2Compute stiffness factor
2
1/ 4 / /
b dT dT
Compute the Kirchoff stress in the fiber fibre i, 2 a
dev
nn
14 Compute the elasticity tensor
, 4 2 / 3
fibre i b dev n n dev nn I I
Update the stress tensor and the elasticity tensor
, fibre fibre fibre i
, ,
fibre i fibre fibre i
ENDIF
For each integration point the collagen resolving rates i and i are computed (note that the disintegration rate is strictly depended of the time t
The undulations min and min is computed and stored in the history vector The collagen density is updated i i
and stored in the history vector ENDO
Note that substantial memory requirements are needs as min , ma x and each
integration point, each gauss point, each element and at each time step t.
5.3 Automatic time stepping
There are several ways of automatically adjusting the size of the time step. One concept is to link the size of the time step to an estimation of the solution error.
5.3.1 Error estimation
In order to verify a solutions correctness, verification and validation is needed. Verification deals with the mathematical correctness of a numerical solution whereas validation deals with the physical correctness of a given model. Consequently, the primary goals of verification is to identify errors in the computer code and to quantify the numerical error in the computed solution. Essentially the theme of verification is error estimation. This thesis will only discuss the verification.
There are three types of numerical errors: round off error, iterative convergence error and discretization error (DE), also known as the truncation error (Roy, 2005). The round off error occurs due to finite arithmetic on digital computers, iterative error occurs due to incomplete convergence and DE is the difference between the exact solution to the discrete equations and the exact solution to the partial differential equations. This thesis will only discuss the DE.
15
,k k
LDE y y t (5.14)
where y is the exact solution at time k t and k y t
k is the approximation of the exact solution at time t . Richardson extrapolation is usually used to estimate the DE. kFor example, for an Forward Euler scheme, the LDE can be approximated by
1
2
1 '' , 2 k k LDE t y t (5.15)which can be shown to be an O
t2 accurate estimation of the error (Kavetski, Binning, & Sloan, 2002).5.3.2 Stability for explicit schemes
Explicit time discretization technique is easy to implement and requires a low cost per time step. An explicit scheme is however unstable and require small time steps.
The stability is not unconditionally. For example, for the Forward Euler method
1 k k, k ,
k
y y hf t y (5.16)
where h is the step size and
k, k
k
dy f t y
dt
,
one can use a model equation according to (5.17) in order to study the linear stability characteristics of the numerical scheme
,
dy y
dt (5.17)
where y
0 y0 and is some characteristic value of the ordinary differential equation (ODE) that arises from assuming that the ODE behaves in a linear manner. We get
1 1 k k k k y y h y y h (5.18)The amplification factor Gk is therefore 1 1 k k k y G h y (5.19)
If we assume that is real, Forward Euler method is stable if 1 1 h 1,
16 which implies
0h2 (5.21)
There is therefore a restriction regarding the size of the time-step h .
The framework presented in 5.1 is non-linear, which implies that is non-linear, and the stability is not given by (5.21). Also, this stability is valid one dimensional problems, the framework in 5.1 is multidimensional. A pragmatic parameter test was conducted to test the automatic time step method presented in 6.2.
5.3.3 Examples of different automatic time stepping schemes
A solution to the time stepping problem is to implement an automatic time stepping technique which changes the value of the time step according to a specified scheme. One example of an adaptive time step scheme is to use a “look-ahead” technique. The outline is presented below Figure 3
Figure 3. Scheme for adaptive time stepping
The LDE can be approximated by (5.15) and should not exceed a specified LDEtolerance . The cost per time step increases substantially, but on the other hand it enhances the robustness and the overall efficiency of the code.
Another simple mechanism is to use an evolutionary proportional-integral derivative (PID) controller where one monitors the relative changes ekaccording to
17 1 1 . k k k k y y e y (5.22)
Here y is an indicator variable, for example in Martuffi and Gassers framework. If ek where tolerance then change the time step using k k
k
t t
e
.
5.4 Basic theory Contact mechanics
An introduction to kinematics and contact constraints is presented in this thesis.
5.4.1 Contact kinematics
Contact occurs either at a point, a line or a surface area. The dimension of the contact surface Γ is one lower than the space dimension. Thus, for a 2D problem Γ is a curve and for a 3D problem Γ is a surface (Kloosterman, 2002). If the bodies, 1
and2 , approach
each other two points 1
X and 2
X in the bodies can occupy the same position in the current configuration
2 1X X , as presented in Figure 4.
Figure 4. Contact illustration, Picture from (Kloosterman, 2002)
To describe the kinematics of the contacting surface a coordinate system on the surface is needed. This is accomplished by a parameterization of the surface, presented in Figure 5.
18
Figure 5. Parameterization of the contact surface, picture from (Kloosterman, 2002)
Two steps are needed to detect if contact takes place or not; the global search and the setup of local kinematic relations.
Enforcing contact constraints yields
2 1
10
x x n (5.24)
where 1
n is the normal vector associated with body 1
and 2
x and x1 are points on the corresponding bodies.
If one assume the contact boundary to be convex, every point 2
x on 2
can be related to a point x1 on 1 by a minimum distance problem
1 1
2 1 2 1 1 2 1 ˆ , min d X x x x x ξ (5.25) The point 1x denotes the minimum distance point on surface 1
in relation to a surface
2
.
Once the minimal distance problem is solved one can define an inequality constraint
2 1
1 Ng x x n (5.26)
which enforces the contact constraints.
For frictionless contact this is the only constraint that needs to be satisfied.
5.4.2 Contact constraint treatments
19
5.5 Discretization of contact constraints
In order to solve complicated contact problems numerical methods such as the Finite Element Analysis (FEA) are used. It is assumed that the reader is familiar with how the method is used to discretize the continuum. This section will describe how FEA is used to discretize the contact constraints (Wriggers, 2006).
5.5.1 LAGRANGE multiplier
For a frictionless case the LAGRANGE multiplier method yields for the residual
1 h c i c n h h N N i N g dN g d
(5.27)where c is the contact boundary, N is the LAGRANGE multipliers, and gN is the gap function. For the constraint equation one obtains
1 0 0 c c h i n h h N N N i N g d g d
(5.28)The discretized interpolations H N and h N g on h ci are defined as
h N K h N I N K I I N K M g N g
(5.29)here is a local convective coordinate. It defines the shape functions NI and MK on the reference element. The nodal value of the LAGRANGE multiplier is NK and gN I is the nodal value of the gap function. The sum in (5.28) is made up until nc active contact constraints. It is not always clear on which side these integrals are supposed to be carried out. This is the reason for defining a master-slave relationship, where one master boundary and one slave boundary is defined.
Using this one obtains a discrete form of the potential energy contributed by contact
LM
according to
,
TLM T T
u A u f A C u (5.30) where u is a displacement vector containing the nodal displacements of both bodies, A is a vector containing all of the LAGRANGE multipliers associated with nc , K is the stiffness matric associated with both bodies, f contains the body forces and surface tractions of both bodies and C is given by
h T
N i
20
5.5.2 Penalty method
In contrast to the LAGRANGE multiplier method which needs discretization in both the field of LAGRANGE multipliers and the displacement field u , the penalty method only needs discretization in the later. The contact constraint is
h c c h N g dN N g dN
(5.32)By performing the integrations the potential energy contribution by contact P
is
2 T P N T u u f CC u (5.33)Variation of virtual displacements yields
, P K K u f (5.34) where P T N K CC .
6. Method
6.1 Geometry and load case
A single brick element was used for testing. Also a representation of an AAA and a representations of an AAA in contact with the spine was considered.
6.1.1 One brick element with prescribed displacement
One Q1-P0 element was prescribed displacement d according to Figure 6. The x-y-plane was fixed in the z-direction.
Figure 6. Geometry and load case 1: Q1-P0 brick
P P
z
z
x
y
21
6.1.2 One brick element with prescribed load
An identical brick described in section 6.1.1 was used, and a force F was prescribed at the same plane as the displacement. The x-y-plane was fixed in the z-direction.
6.1.3 Abdominal Aorta with internal pressure
A representation of an AAA, see Figure 7, consiting of 2 503 Q1P0 elements and 2561 nodal points. An internal pressure of 13.3 kPa was prescriped on the internal surface. The model was fixed at the bottom and top slices (nodal degrees of freedom were locked).
The internal pressure was proportionally applied during n iterations according to Figure 8 and thereafter held constant until the end of the simulation. This was required due to the non-linear behavior of the model.
6.1.4 Abdominal Aorta with internal pressure and contact with spine
A digital representation of the geomentry was modeled from CT-scans shown in Figure 9.
The AAA was modeled using two materials: the vessel wall conisting of Q1P0 elements, using the constitutuve framwork presented in 5.1, and a intra-luminal thrombos, ILT. A pressure of 13.3 kPa was prescriped on the internal surface. The model was fixed at top slices (nodal degrees of freedom were locked) and free in the axial direction on the bottom slice.
The spine was modeled using quadratic elements and it was fixed in all degrees (rigid). The digital representation is presented in Figure 10.
The contact was modelded using the penalty method with 10
(6.1)
The spine was choosen to be the slave surface and the entire surface of the AAA vessel wall the master surface.
22
Figure 7. Geometry AAA
Figure 8: Proportional load, load case AAA with internal pressure. The load 13.6 kPa was constant until end of simulation. 0 2 4 6 8 10 12 14 1 2 3 4 5 6 7 8 9 n end kPa Iteration
23
Figure 9. CT scan of the AAA and spine
Figure 10. AAA and spine, von Mises stress.
Trombos
Wall
24
6.2 Automatic time stepping
In this project a technique which connects the size of the time step with the collagen dynamics. The idea behind this is that the collagen production rate is implicitly connected to the accuracy of the solution.By examining equation(5.5) and section 5.1.2 one can notice that
..., min, max
...,
fibre f f
(6.2)
Two different approaches were tested. The first in which the largest incremental change in density in of directions was considered according to
1 2 , ... n MAX (6.3)
where the sub-indices denotes the spatial direction of the collagen fiber was aligned. The approach summed the incremental change in density in all directions according to
1 , n i i
(6.4)and the incremental change in min and ma x also was taken in to account.
25
Figure 11. Outline of the automatic time stepping scheme
The solution command to activate the automatic time stepping is
1
2
3AUTO MATErial rvalu rvalu rvalu (6.5)
where rvalu(1) – (3) are parameters used to control the value of rmeas. The value rmeas is passed between the constitutive model and the solution by a common statement according to
𝑟𝑒𝑎𝑙 ∗ 8 𝑟𝑚𝑒𝑎𝑠, 𝑟𝑣𝑎𝑙𝑢 𝑙𝑜𝑔𝑖𝑐𝑎𝑙 𝑎𝑟𝑎𝑡𝑓𝑙
(𝑐𝑜𝑚𝑚𝑜𝑛/𝑒𝑙𝑎𝑢𝑡𝑜)/𝑟𝑚𝑒𝑎𝑠, 𝑟𝑣𝑎𝑙𝑢(3), 𝑎𝑟𝑎𝑓𝑙
6.2.1 Automatic time stepping first approach
In the first approach the value rmeas is connected to the collagen production rate by
max , , rmeas rmeas rvalu (6.6)
where rvalu is a user specified value used to control rmeas. The specified value rvalu essentially governs the maximum accepted value of rmeas.
26
6.2.2 Automatic time stepping second approach
In the second approach the value rmeas in connected to both the net contribution of collagen, in all directions, and the change in in min and ma x according to
max max( , ) / rvalu(1), ),
rmeas rmeas (6.7) where 1 / , t t tot tot (6.8) and 1 1
min min max max
100 t t t t
(6.9)
The time step is thus dependent on two structural parameters.
6.2.3 Estimation of error
The error estimation was carried out by a parameter test according to Figure 12. The test was carried out both for the geometry and load case described in 6.1.1 and 6.1.3. The test was carried out up to a specified time t and a specified constantt.
Figure 12. Outline of the parameter test
In order to compare the effectiveness of the two time-stepping methods the elapsed CPU-time was obtained by the command line
27
where CPU_TIME returns a real value of the elapsed CPU-time in seconds. The elapsed CPU time is offset by the system clock which renders an absolute CPU-time useless. The offset is written in the output file which FEAP provides. The comparison value of CPU time is obtained by decreasing the absolute CPU time with the offset provided by FEAP according to (6.11)
time abs offset
CPU CPU CPU (6.11)
The error is estimated by
i,ma i,max i,max, x, , , r e f i er r or ref (6.12)
where i m a x, is the maximum principal Cauchy stress calculated by the automatic
time-stepping in element i and i,max,r e f is the maximum Cauchy stress calculated by a
constant t in element i. The element i is an element at the location of large stress. The Cauchy stress were evaluated at the Gauss point level and averaged over the element. One should note that the reference stress i,max,r e f is not necessarily the exact solution
but regarded as a good estimation of the exact solution.
7. Results
Results are presented below. The first time stepping approach described in 6.2.1 was used for the two initial sections 7.1.1 and 7.1.2. In contrast the following sections, 7.2.1 and 7.2.2, used the second approach presented in 6.2.2.
7.1 First time stepping approach
28
Figure 13: Parameters for the first automatic time stepping scheme
7.1.1 Results: One brick element with prescribed displacement
The simulations with the one brick element were carried out for 5 different t (from 3,25 days to two months of real physical time), see Table 2. The pseudo time is the sum of the time steps used in FEAP.
Table 2: Computational and physical time, load case 1
Pseudo time t [1] Physical time (days)
1 3.25
2 7.5
4 15
8 30
16 60
The resulting first principal stress and CPU-time, using the specified constant time step listed in Figure 13, are presented in Figure 14. The stress relaxes due to optimization of the collagen structure.
29
Figure 14. Maximum principal stress and CPU-time, load case 1, constant time step
Two different R-values were used for the automatic time stepping, according to Figure 13.
The resulting relative error relative to the constant time step t 0.001 is presented in Figure 15.
Figure 15. Relative error for different R-values, load case 1, compared to a constant time step of 0.001.
The CPU-time for the automatic time-step scheme is presented in Figure 16. 0 50 100 150 200 250 300 350 400 450 0.00E+00 0.10E-03 0.20E-03 0.30E-03 0.40E-03 0.50E-03 0.60E-03 0.70E-03 0.80E-03 0.90E-03 1.00E-03 3,25 7,5 15 30 60
Phycial time [days]
CPU -time [ s] Principa l s tre ss [ MPa]
Maximum princial stress and CPU-time for constant time
step
Maximum principal stress [MPa] CPU-time [s]
0 0,5 1 1,5 2 2,5 3,25 7,5 15 30 60 Re lativ e Error [% ]
Physical time [days]
Relative error, principal Cauchy stress 1, for different
R-values
R-value = 1.0E-3 R-value =1.25E-3
30
Figure 16. CPU-time for different R-values, load case 1.
Comparing at t16 (60 physical days) the automatic time step scheme is much faster, with a small relative error, see Table 3.
Table 3: Results at t=16 (60 physical days), automatic and constant time step
Constant 0.001 t
R-value 1E-3 R-value 1.25E-3
Percent speed-up - 19 967% 18 350%
Relative error - 0.21-0.43% 0.21-0.43%
7.1.2 Results: Abdominal Aorta with internal pressure
The results were compared at element 802, which lies in the area of large stress, see Figure 17. 0.0 0.5 1.0 1.5 2.0 2.5 3,25 7,5 15 30 60 C PU -tim e [s ]
Physical time [days]
CPU-time for different R-values
31
Figure 17: Load case 2, results, element 802
First, a constant t was chosen to 0.001. This was the largest t to converge at time t. The simulations were run up to the pseudo timet 1 , which corresponds to 3.25 days of growth. The computational load made it difficult to run the test to a larger t due to time restraints.
The average principal Cauchy stress-values is presented in Table 4. The computation took 10 424 CPU-time.
Table 4: Stress and CPU-time, load case 2, constant time step, t 0.001
Maximum principal stress [MPa]
Principal stress 1 0.168
Principal stress 2 0.0155
Principal stress 3 -7.38E-3
For the automatic time-stepping scheme the R-values tested were R=1.0E-5 and 1.25E-5. Larger R-values did not converge. Also, t values larger or equal than 1.81 did not converge due to limitations in the collagen rate production (see discussion).
32
The relative error compared to the constant time step of t 0.001 and pseudo time 1
t (3.25 days) for the two different R-values are presented in Figure 18.
The automatic time step scheme was significantly faster compared to a constant time step, with a marginal relative error, see Table 5.
Table 5: Comparison between constant time step and automatic at t=1 (3.25 days)
Constant t 0.001 R-value 1E-5 R-value 1.25E-5
Percent speed-up - 50% 74%
Relative error - 0.97-1.1% 0.97-1.1%
7.2 Second time stepping approach
The maximum time step was set to tma x 100 for both load cases.
7.2.1 Results: One brick element with prescribed load
The displacement for each time step is presented in Figure 19. The initial collagen structure is not optimized which yields the curve seen up to pseudo time = 150. Afterwards the element reaches homeostasis.
5400 5600 5800 6000 6200 6400 6600 6800 7000 7200 0,85% 0,90% 0,95% 1,00% 1,05% 1,10% 1,15% R=1E-5 R=1.25E-5 CPU -time [ s] Re lativ e error
Relalative error compared to constant Δt=0.001 and CPU
time for two different R-values
Relative error: Principal stress 1 Relative error: Principal stress 2 CPU-time [s]
Figure 18. Relative error and CPU-time, AAA with internal pressure.
33
Figure 19. Time vs displacement, one brick element with prescribed load.
The size of the time-step for each time-step is presented in Figure 20. The characteristic erratic behavior of the incremental size is due to the time step being dependent on two different structural parameters.
Figure 20. Time vs size of time-step, one brick element with prescribed load.
7.2.2 Results: Abdominal aorta with internal pressure and contact with spine
Two different values for the maximum collagen production max
34 The size of the time-step vs time for max 1.2
and max 0.7
is presented in Figure 21 and Figure 22. The value of the maximum collagen production significantly influenced the growth of the AAA.
Using max 1.2
the AAA was stable and did not reach failure. Using max 0.7
the AAA reached failure at t=300 (1125 physical days).
Figure 21. Size of time step vs time, AAA and spine, rhomax = 1.2.
Figure 22. Size of time-step vs time, AAA and spine, rhomax = 0.7.
35
Figure 23. Penetration of spine, showing von Mises stress [kPA].
For both cases, the AAA experienced a stress maxima on one side, and middle of the AAA. The area of low stress in connection to the stress maxima is due to the ILT carries load, see Figure 24.
Figure 24. Von Mises stress [kPa], AAA and spine, min and max, rhomax = 1.2.
8. Discussion and Conclusions
The objective of the thesis was to implement an automatic time stepping scheme and test its validity. Also, the effect for the stress distribution in the AAA resulting from contact with the spine was supposed to be tested.
The automatic time-step scheme showed significantly better computational performance compared to a fixed time step. This was achieved with a small relative error. One should however note that the reference stress, computed using a fixed time step, i,max,r e f is not necessarily the exact solution but regarded as a good estimation of
the exact solution.
36
The thesis presents a way of connecting the size of the time step to the rate of the collagen production and how the collagen is added to the existing distribution. This enables the user to control the performance/error ratio of the calculation due to the proportionality between the collagen rate production and the correctness of the solution.
The reason for non-convergence for larger or equal values of t=1.81, for the first automatic time-stepping scheme, was due to limitations of the model of collagen production (Gasser C. , 2014). The size of the time step using the suggested scheme is essentially inversely proportional to the size of the rate of the collagen production, see (6.6). This means that if the collagen production rate grows the time step gets smaller. This is exactly what the program experienced – the time step approached zero for large values of t.
The second automatic time-stepping scheme described in 6.2.2 is dependent on two different structural parameters. This is probably the reason for the oscillating behavior of the size of the time step seen in Figure 20 and Figure 21.
The appearance in Figure 22 is also expected. The collagen structure was not optimized initially. The AAA experienced shrinking the first time-steps and small incremental change in net collagen addition, thus the time step was allowed to increase. When the AAA starts expanding, more collagen is needed to bear the load. A large contribution of collagen per time step decreases the size of the time-step. The size of the time-step decrease until failure is reached. Interesting enough, the stress at which failure occurred was in the same order of magnitude as when rupture usually occurs (Vorp & Vande Geest, 2005).
In order to verify the correctness of the numerical scheme significant work is still do be done. Especially for larger values of t. The second automatic time-scheme should be validates in the same manner as the first. Also, in order to validate the model, the results should be compared to scans of in vitro diseased aortas at set intervals.
Other automatic time stepping schemes could also be tested. For example a “look ahead” technique. Even though the cost per time step I substantially higher, significant increase in performance could possibly be obtained.
Also, the variables controlling when to increase or decrease the time-step, for given values of RMEAS, could be varied in order to optimize the scheme for the geometry of the aorta and the constitutive model.
Lastly, optimizing the penalty value in the contact simulation, should be investigated and optimizing the contact surfaces (restricting them to the area of contact).
9. References
37
Gasser, T., & Martufi, G. (2011). A Constitutive model for Vascular tissue that intergrates fibtil, fiber and continuum levels with application to the isotropic and passive properties of the infrarenal aorta. Journal of Biomechanics, 2544-2550.
J.Burkart. (2014, 1 28). SPHERE_DESIGN_RULE. Retrieved from Hardin and Sloane
Spherical Designs:
http://people.sc.fsu.edu/~jburkardt%20/f_src/sphere_design_rule/sphere_desi gn_rule.html
Kavetski, D., Binning, P., & Sloan, S. W. (2002). Adaptive backward Euler time stepping with truncation error control for numerical modelling of unsaturated 3uid 3ow. INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, 1301-1322.
Kloosterman, G. (2002). Contact Methods in FInite Elment Simulations. Haag: Netherlands Institute for Metals Research.
Martufi, G., & Gasser, C. T. (2012). Turnover of fibrillar collagen in soft biological tissue with application to the expansion of abdominal aortic aneurysms. J. R. Soc. Interface , 3366-3377.
Roy, C. J. (2005). Review of code and solution verification procedures for computational simulation. Journal of Computational Physics, 131-156.
Rupp, C., Howard, M., & Weickum, G. (2006). Incompressible Mixed (u/p) Elements for the CAS FEM Cod. Boulder: Center for Aerospace Structures Department of Aerospace Engineering Sciences University of Colorado at Boulder.
Simo, C. J., & Taylor, L. R. (1985). Consistent tangent operators for rate-independent elastoplasticity'. Computer Methods in Applied Mechanics and Engineering, 101-118.
Taylor, L. R. (2011). FEAP -- A Finite Elemement Analysis Program, 8.3 User Manual. Berkeley: Department of Civil and Enviromental Engineering University of California at Berkley.
Taylor, L. R. (2013). FEAP - A Finite Element Analysis Program Version 8.4 Programmer Manual. Berkley: University of California at Berkley.
Upchurch, G. R., & Schaub, T. A. (2006). Abdominal Aortic Aneurysm. Am Fam Physician, 1198-1204.
38
Vorp, D. A., & Vande Geest, P. J. (2005). Biomechanical Determinants of Abdominal Aortic Aneurysm Rupture. Arteriosclerosis, Thrombosis, and Vascular Biology., 1558-1566.
