Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

52, 4, pp. 917-925, Warsaw 2014

NUMERICAL SIMULATIONS OF ARTERIES WITH AN ADAPTIVE

FINITE ELEMENT METHOD

A. Karolina Fuksa
Cracow University of Technology, Departament of Civil Engineering, Kraków, Poland

e-mail: afuksa@pk.edu.pl

Waldemar Rachowicz
Cracow University of Technology, Institute of Computer Science, Kraków, Poland

e-mail: wrachowicz@pk.edu.pl

This paper presents applications of the adaptive FEM in computational biomechanics. The
study involves finite deformation of artery walls built of layers of elastic, nearly incompres-
sible (rubber-like) materials. The models are reinforced with nearly inextensible collagen
fibres stretched along the directions close to the circumferential direction. A Total Lagran-
gian technique based on the displacement-pressure formulation is presented. Residual error
estimates and finite element mesh adaptation strategies are developed. The procedures are
illustrated for solutions of arteries under various load considerations.

Keywords: computational biomechanics, finite elementmethod, adaptivity, error estimation,
finite deformation

1. Introduction

Computational biomechanics is a rapidly developing branch of numerical simulations based on
the Finite ElementMethod approximation. Research in themechanics of soft tissues and parti-
cularly problems related to the simulations of aortic aneurysms and arteries with changes due to
atherosclerosis are of particular interest. In both cases, the arteries are subjected to a high risk
of rupturewhenpressurized.A commonly accepted approach for such biomechanical simulations
is the use of low-order, very dense finite element (FE)meshes to discretize the biological objects
of interest. Such an approach ismotivated by the necessity for accurate approximation of a com-
plex, patient-specific geometry involving a part of the artery. Themain difficulties in obtaining
a reliable simulation often lie in the ability to determine an appropriate constitutive relation
between the stress and strain and also in the assessment of the material parameters through
measurements. Representing accurately the mechanical behaviour of the real-life tissue is thus
considered the main difficulty. The approximation error related to the discretization errors is
seldom estimated or suppressed by successive mesh refinements, i.e. by convergence studies.
However, there exists a well established school of computational mechanics which showed

that the discretization error can be controlled without loss of effciency (Demkowicz et al., 2007).
One can use relatively coarse FE meshes and adapt them using h-refinements, i.e. changing
locally the size of elements h, and/or p-refinements which involves changing locally the order
of approximation p. The mesh adaptations are intended to reduce the estimated discretization
error below a prescribed threshold for a minimum computational effort. The estimated error
might be understood globally or it may concern some user-specified quantities of interest such
as maximum equivalent stress. Such procedures are most often used for problems with simple
andwell-knownphysicalmodels of thematerial. Thequestion then ariseswhether the techniques
of error estimation andmesh adaptivity could represent an advantage in computational biome-
chanics. Although themodeling error due to poormaterial descriptionmay be considered as the
main cause of poormechanical response predictions today, the situationmay change so that the



918 A.K. Fuksa,W. Rachowicz

discretization error becomes important. Firstly, new effective versions of the finite element me-
thod arise such as the iso-geometric analysis (IGA) or theDiscontinous Petrov-Galerkin (DPG)
techniqueswhichallow anaccurate approximation of solutions on complex geometries using rela-
tively coarse meshes. Secondly, the constant progress in the development of constitutive models
and the recovery of their parameters may reduce in the future the dominance of this factor as a
source of inaccuracy for simulations. Also, the interaction of biological tissueswithmedical devi-
ces such as stentsmay require error estimation. Finally, higher order approximation is in general
more resistant to numerical locking which might be caused by the near-incompressibility and
near-inextensibility of the arterial tissue reinforced with collagen fibres (Holzapfel and Gasser,
2001).

2. Adaptive FEM for the mechanics of arteries

In this study, we apply the methodology of error estimation and adaptivity for finite element
simulations of arterial mechanics. Appropriate algorithms have already been developed and
some adjustments have been made to obtain an accurate simulation. We developed the two-
-field formulation of displacement-pressure type (Simo and Taylor, 1991) with no inter-element
continuity enforced for the pressure variable. The Total and Updated Lagrangian approaches
can be used to solve the finite elasticity problem. We could also use an automatic h-adaptive
algorithm (Demkowicz et al., 2007) to generate adaptive meshes based on the adjustment of
the mesh to reduce the interpolation error of the solution. We developed an error-estimation
technique which is applicable to finite elasticity problems. It is a version of the residualmethod
proposed by Bank andWeiser (1985) in the context of linear elliptic boundary value problems.
It may serve for error estimation of the solution in a global (energy) norm as well as errors of
local quantities of interest in its goal-oriented version. The distribution of error indicators may
also be used to guide adaptivity of finite element meshes. The adaptive FEM has been applied
to solve typical problems in arterial mechanics.
A typical artery consists of concentric layers of tissue forming a tubular structure with

properties significantly varying across the thickness. Such composites are nearly-incompressible
and usually described as neo-Hookean materials reinforced with two families of collagen fibres.
Atherosclerosis is a vascular disease which causes narrowing of the artery internal diameter
and hardening of its walls as a result of plaque deposits. In some treatment therapies, a stent
is inserted into the artery and is adequately deformed to widen the cross-section narrowed by
atherosclerosis. Aprecise control of this treatment called angioplasty is necessary. It is a delicate
procedurewhich requires sufficient widening of the arterywithout causing damage to the artery
walls (Edelman andRogers, 1998). In this research, we perform preliminary tests on amodel of
an artery under various load considerations for further investigations on the interaction between
medical devices and biological tissues like for the case of angioplasty. Simulations of this sort
involve a contact problemwith large sliding andwith finite deformation of the bodies in contact.
The reliability of numerical simulations is also essential following the growth and remodelling
(G&R) of aneurysms. To start, we wish to assess the danger of rupture of an in vitro artery
under variousmechanical loads. Investigations of this kind involvemodeling of a patient-specific
geometry of theartery forwhich IGA is an interesting option. Such researchbecomes increasingly
important for medical applications and the treatment of cardiovascular conditions.

3. Rubber-like materials reinforced with fibres

The tissue of a typical arterymaybe considered as a hyperelastic nearly-incompressiblematerial
reinforcedwith two families of collagen fibres (Holzapfel andGasser, 2001) forwhich the stiffness
grows exponentially with stretch in a way that they become nearly inextensible. The fibres are



Numerical simulations of arteries with an adaptive finite element method 919

inactive if theyare subject to compression.Acommonway to characterize ahyperelasticmaterial
consists of postulating its internal strain energyas a functionof the invariants of thedeformation.
In thematerial description of the deformation of a body x=x(X, t), X and x are the locations
of amaterial particle in the initial and deformed configurations respectively, we define the right
Cauchy-Green deformation tensor

C=FTF (3.1)

where F = ∂x/∂X is the deformation gradient. Defining by J = det(C)1/2 = det(F) > 0 the
volume ratio, we can introduce the multiplicative decomposition of the deformation measures
F and C

F= J1/3F C= J2/3C (3.2)

where F and C are the unimodular deformation gradient and the Cauchy-Green tensor, re-
spectively. We assume two families of fibre directions ai(X), |a|= 1, i= 1,2 which define the
structural tensors

Ai =ai⊗ai i=1,2 (3.3)

With the deformation measures defined above, we associate the following set of invariants

I1 = tr(C) I2 =
1

2
[(trC)2− trC2] I4i =C :Ai i=1,2 (3.4)

and their unimodular counterparts

I1 = tr(C) I2 =
1

2
[(trC)2− trC

2
] I4i =C :Ai i=1,2 (3.5)

Further, we assume an additive split of the internal strain energy function into the volumetric,
isochoric and fibre components

W =Wvol(J)+Wiso(I1,I2)+
2∑

i=1

Wfi(I4i) (3.6)

Splitting the internal energy into volumetric and isochoric parts for nearly incompressible iso-
tropic materials is attributed to Flory (1961). The constitutive relation for the second Piola-
Kirchhoff stress tensor S, which is implied by (3.5), reads

S=−pJC−1+2J−2/3
{
Wiso,1Dev[I]−Wiso,2Dev[C

−2
]+

2∑

i=1

Wfi,4Dev[Ai]
}

(3.7)

where p denotes the hydrostatic pressure

p=−
∂Wvol
∂J

(3.8)

and Dev[•] = [•]−1/3(C : (•))C−1. We assume the following form of the subsequent compo-
nents of W (Holzapfel and Gasser, 2001)

Wvol =
K

4
(J−1)2 Wiso =

µ

2
(I1−3) Wfi =

k1j
2k2i

{
exp[k2i(I4i−1)

2]−1
}
(3.9)

where K, µ, k1i, k2i represent the parameters characterizing the material.



920 A.K. Fuksa,W. Rachowicz

4. Formulation of the problem and finite element approximation

Typically, the parameter K is significantly larger than µ, reflecting the fact that the material
is nearly incompressible. It becomes convenient to formulate a two-field mixed problem with
displacements, u= x−X, and the pressure p as the unknowns. This involves the principle of
virtual work and the weak form of the constitutive formula for the pressure (Simo and Taylor,
1991)

∫

Ω

S :
1

2
(FT∇Xv+∇

T
XvF) dX =

∫

Ω

Jb ·v dX+

∫

ΓN

(FTv) · t̂ dA ∀v∈V

∫

Ω
−
( 2
K
p− (J−1)

)
q dX =0 ∀q∈L2(Ω)

(4.1)

In this statement, b denotes the forces per unit volume, t̂ are the tractions prescribed
on the Neumann part of the boundary ΓN, v and q are appropriate test functions with
V = {v ∈ H1(Ω) : v = 0 on ΓD} and ΓD stands for the part of ∂Ω with Dirichlet bounda-
ry conditions. This nonlinear problem can be linearized and solved with the Newton-Raphson
procedure. The linearized problem takes the following form

a(∆u,v)+ b(∆p,v)= l1(v) ∀v∈H
1(Ω)

b(q,∆u)− c(∆p,q)= l2(q) ∀q∈L
2(Ω)

(4.2)

where the bilinear forms are defined as follows

a(∆u,v) =

∫

Ω

E
′(u,v) : [Cdev−pJ(C

−1⊗C−1−2C−1⊙C−1)] :E′

+∇Xv : [∇X∆u(Dev[S]−pJC
−1)]dX

b(∆p,v) =−

∫

Ω

J∇X ·v∆pdX

(4.3)

where

∇X ·v=C
−1 :E′(u,v) E′(u,v) =

1

2
(∇TXvF+F

T
∇Xv) (4.4)

and

Cdev =2
∂

∂C

{
J−2/3Dev

[
2
∂W

∂C

]}
(4.5)

denotes the tensor of elasticities. The linear functionals li(v) correspond to the residuumof the
current approximation (u,p). To discretize the problem hexahedral elements of the Qp/Pp−1

family are used where p is the order of approximation. It provides a stable approximation of
the mixed problem.

5. Error estimation and adaptivity

Estimation of the approximation error of themixed formulation of finite elasticity problem (4.1)
was considered by Rüter and Stein (2000). It is based on the Helmholtz decomposition of the
error in the first Piola-Kirchhoff stress tensor, P−Ph =∇ψ+∇×φ, where ψ and φ are the



Numerical simulations of arteries with an adaptive finite element method 921

error functions. Following this representation, the residual of equilibrium equation (4.1)1 can be
written as follows

R(v) =

∫

Ω

Jb ·v dX+

∫

ΓN

t̂ · (FTv) dA−

∫

Ω

Ph :∇v dX

=

∫

Ω

(P−Ph) :∇v dX =(ψ,v)1

(5.1)

which defines the function ψ as a solution to the variational boundary-value problem

Find ψ : (ψ,v)1 =R(v) ∀v∈V (5.2)

Thenormof ψ is estimatedusing the technique of error estimationwith self-equilibrated element
residuals due to Ainsworth and Oden (1992). The final estimate of the error in displacements
and in pressure is found as follows

|u−uh|1,Ω +‖p−ph‖0,Ω ¬C
(
|ψ|1,Ω +

∥∥∥U −
2

K
ph
∥∥∥
0,Ω

)
(5.3)

where U =J−1. That is, it includes the norm of the error function ψ and themeasure of the
dissatisfaction of the constitutive relation for pressure.
Wemodify this approachbyreplacing themethodof self-equilibrated residualsby theelement

residual technique proposed by Bank and Weiser (1985). In this technique, we assume that
Vhp(K) is the space of shape functions of the element K. We consider an enriched version of
this space Vh,p+1(K), which is obtained by enriching the order of approximation by one. Let
Πhp be an interpolation operator defined for Vhp(K), for instance of the Lagrangian type. We
consider a kernel of this operator in Vh,p+1(K)

MK = {v∈Vh,p+1 : Πhpv=0} (5.4)

For every element K, we solve a discrete local boundary-value problem for the error indicator
function φK

φK ∈MK : a(φK,v) =RK(v) ∀v∈MK (5.5)

where RK are the element residuals and a(·, ·) is a bilinear form of the problem. The estimate
of the dual norm of the global residual is expressed as follows

‖R‖∗ ¬C

(
∑

K

‖φK‖
2
E,K

)1/2
(5.6)

where ‖ · ‖E is the energy norm associated with the form a. An advantage of this technique is
that it is reliable and inexpensive as the size of the local problems is O(p2) instead of O(p3) for
other techniques including themethod of self-equilibrated residuals.
The distribution of error indicators is a basis for local refinements of the finite elementmesh.

A standard feedbackmesh adaptation algorithmdivides the elements with the largest errors and
finds the solution on the newmesh, and so on. It can be stated as follows:

1. Solve the problem on the current mesh.

2. Find the error indicators ηK.

3. Stop if (
∑
K η
2
K)
1/2 ¬TOL.

4. Refine the elements K for which: ηK >αmaxLηL.

5. Go to 1.

Here, 0<α< 1 is a user-defined parameter controlling the adaptive process.



922 A.K. Fuksa,W. Rachowicz

6. Numerical examples

6.1. Yosibash and Priel benchmark

In the first example, we consider the effect of pressurizing a typical artery such as the left
anterior descending coronary artery (LAD)which consists of twomain layers, the internalmedia
and the external adventitia layer.Thisproblemwas investigated as abenchmarktest byYosibash
and Priel (2011). The geometry of the artery can bemodeled as a cylindrical tube composed of
two layers of the material. The material data for the media is K = 100000kPa, µ= 54.0kPa,
k1 =0.64kPa and k2 =3.54 and the fibres are inclined with respect to the horizontal direction
at an angle of ±10◦. In the adventitia layer, K = 100000kPa, µ = 5.4kPa, k1 = 5.1kPa
and k2 = 15.4 and the fibres are inclined with respect to the horizontal direction at an angle
of ±40◦. The artery is pressurized by an internal pressure of 13.3kPa (100mmHg). Tetrahedral
elements with the Q2/P1 shape functions are used. Figure 1 shows the geometry of the artery,
the deformed configuration as well as the pressure distribution.

Fig. 1. Yosibash and Priel benchmark: (a) classes of the anisotropic material, (b) radial displacements,
(c) pressure along the layers of media and adventitia

6.2. LAD aneurysm

In the following test, amodel of the artery similar to the one used in theYosibash andPriel
benchmark is considered. However, an additional cavity within the media layer is created to
model themedical condition of aneurysm.Thematerial data remains the sameas in the previous
example and the artery is subjected to an internal pressure of 13.3kPa. The distribution of the
deviatoric stress and pressure on the surface of the artery are shown in Fig. 2.

6.3. Carotid bifurcation

A typical geometry of an artery such as a carotid bifurcation is depicted in Fig. 3. Also
shown are the classes of thematerial, themedia and the adventitia layers for which thematerial
parameters are K =100000kPa,µ=35.74kPa,k1 =13.9kPa,k2 =13.2.Thefibres are inclined
with respect to the circumferential direction at an angle of ±17.8◦ in the media layer and an
angle of ±30.1◦ in the adventitia layer. In addition, the directions of the fibres are orthogonal in
the area of bifurcation and parallel in the patch located in the vicinity of the saddle point. The
artery is pressurized with an internal pressure of 5.5kPa. Figures 3b and 3c show the pressure
distribution over the external and internal surface of the artery model. The distribution of the
deviatoric stress over the external surface is shown in Fig. 4.



Numerical simulations of arteries with an adaptive finite element method 923

Fig. 2. LAD aneurysm: (a) Dev[S]xx, (b) Dev[S]yy, (c) pressure along the external surface of the
adventitia

Fig. 3. A carotid bifurcation under internal pressure: (a) classes of the anisotropic material,
(b) distribution of pressure on the external and (c) internal surface of the artery

Fig. 4. A carotid bifurcation under internal pressure: (a) Dev[S]xx, (b) Dev[S]yy, (c) Dev[S]zz



924 A.K. Fuksa,W. Rachowicz

6.4. Contact with a rigid body

The future goal of this research is the investigation of the interaction between an artery and
amedical device. It involves a contact problem of deforming bodies. In this Section, we present
a preliminary test of contact with a rigid obstacle.
In this problem, the contact of an elastic cylinder of radius r = 8cm and a rigid obstacle

such as a cylinder of radius r = 20cm is considered. The material parameters of the cylinder
are K =100000kPa,µ=1.0kPawithout fibres in this case. Only normal contact is considered.
The normal gap function is defined as follows. Given a point of the elastic cylinder located at
x2 we find its closest point x1 on the obstacle and the unit outward normal on this surface n1.
The gap function is set as

gN =(x
2−x1) ·n1 (6.1)

No penetration of the bodies is expressed as gN ­ 0, with an additional condition that the
normal stress on the contact area is negative if gN =0,while for gN > 0 it vanishes.We enforce
these conditions by adding to the load side of formulation (4.1)1 the penalty term

∫

Γc

1

ǫ
gNδgN dA (6.2)

where ǫ is a small penalty parameter and Γc is the contact area. The penalty term is activated
only if gN ¬ 0. The variation δgN =v ·n

1 where v stands for a variation of displacements (a
test function). Appropriate linearization of the penalty term is necessary in order to solve the
nonlinear problemof finite elasticity with contactwith theNewton-Raphson iteration (Wriggers,
2002). Selected characteristics of the solutionarepresented inFig. 5.Asimilarproblemof contact
of an elastic sphere with a horizontal rigid plane is illustrated in Fig. 6.

Fig. 5. Contact of the cylinder and a rigid obstacle: (a) displacements uz, (b) Dev[S]zz

Fig. 6. Contact of an elastic sphere with a rigid plane: (a) displacements uz, (b) Dev[S]zz



Numerical simulations of arteries with an adaptive finite element method 925

7. Conclusions

This study investigates the use of the adaptive version of the FEM in the biomechanics of ar-
teries. In particular, the adaptive FEM has been applied to simulate the behaviour of a typical
artery subjected to internal pressure, the pressurization of an artery with aneurysm as well as
a contact problem. Such simulations serve as preliminary studies for further research on the
interaction between amedical device and biological tissues. Thematerial of the artery has been
modeled as nearly-incompressible, neo-Hookean composite reinforced with two families of colla-
gen fibres. TheTotal Lagrangian technique based on the displacement-pressure formulation has
been presented as well as residual error estimate and finite element mesh adaptation strategies.
Finally, an adaptive FEM has been applied to solve typical problems in arterial biomechanics.

Acknowledgement

This work was supported by the Polish National Science Centre under grant number UMO-

2011/01/B/ST6/07306

References

1. Ainsworth M., Oden J.T., 1992, A procedure for a posteriori error estimation for h-p finite
element methods,Computer Methods in Applied Mechanics and Engineering, 101, 73-96

2. Bank R.E., Weiser A., 1985, Some a posteriori error estimates for elliptic partial differential
equations,Mathematics of Computation, 44, 283-301

3. DemkowiczL.,Kurtz J., PardoD.,PaszyńskiM.,RachowiczW.,ZdunekZ., 2007,Com-
puting with hp Finite Elements. II. Frontiers: Three-Dimensional Elliptic and Maxwell Problems

with Applications, CRCPress, Taylor and Francis

4. Edelman E.R., Rogers C.R., 1998, Pathobiologic responses to stenting, American Journal of
Cardiology, 81, 4E-6E

5. FloryR., 1961,Thermodynamic relations for highly elasticmaterials,Transactions of the Faraday
Society, 57, 829-838

6. Holzapfel G.A., Gasser Th.C., 2001, A viscoelastic model for fibre-reinforced composites at
finite strains: Continuum basis, computational aspects and applications, Computer Methods in
Applied Mechanics and Engineering, 190, 4379-4403

7. Rüter M., Stein E., 2000, Analysis, finite element computation and error estimation in trans-
versely isotropic nearly incompressible finite elasticity, Computer Methods in Applied Mechanics
and Engineering, 190, 519-541

8. Simo J.C., TaylorR.L., 1991,Quasi-incompressible finite elasticity in principal stretches,Com-
puter Methods in Applied Mechanics and Engineering, 85, 273-310

9. Wriggers P., 2002,Computational Contact Mechanics, JohnWiley & Sons, 520p.

10. YosibashZ., PrielE., 2011, p-FEMs for hyperelasic nearly-incompressiblematerials under finite
deformationswith applications to arteries simulation, International Journal for NumericalMethods
in Engineering, 88, 1152-1174

Manuscript received October 19, 2013; accepted for print April 11, 2014