Arterial Wall Mechanics

Arteries are blood vessels that carry freshly oxygenated blood from the heart throughout the rest of the body. They are layered structures with the intima inside, followed by the media, and the adventitia. The two outer layers are predominantly responsible for the mechanical behavior of healthy arteries. Both layers are made of collagenous soft tissues that show prominent strain stiffening. Families of collagen fibers give each layer anisotropic properties. These fiber reinforced structures enable blood vessels to sustain large elastic deformation.

The Holzapfel-Gasser-Ogden (HGO) constitutive model described in Ref. 1 captures the anisotropic nonlinear mechanical response observed in experiments on excised arteries.

This model demonstrates how this hyperelastic material is implemented in COMSOL Multiphysics, and the results are compared to those reported in Ref. 1.

Model Definition

The model geometry represents a sector of a carotid artery from a rabbit. Following Ref. 1, the media and adventitia are modeled as a layered cylindrical tube. Model symmetry allows the use of 2D axisymmetric model. Main dimensions are reported in Figure 1.

Figure 1: Carotid artery section of L = 2.5 mm in length. The inner radius Ri is 0.71 mm, the outer radius Ro is 1.1 mm, the media thickness is 0.26 mm, and the adventitia thickness is 0.13 mm.

Typical mechanical experiments measure the response of arterial sections subject to combined axial stretch and internal blood pressure. The set of boundary conditions replicate these experiments.

The symmetry boundary condition allows the bottom end of the artery to freely expand in the radial direction. On the top surface, prescribed displacements in the axial direction account for the axial stretching. Internal pressure is applied with a pressure boundary load on the inner surface.

This model considers axial stretches between 1.5 and 1.9, and internal pressures between 0 and 160 mmHg. The mechanical response in this range is highly nonlinear, resulting in large elastic deformations, and it is mathematically described within the theory of hyperelasticity.

The HGO model is an incompressible anisotropic hyperelastic material model defined by an isochoric strain energy density.

The incompressibility condition implies adding a volumetric stress Svol.

The auxiliary pressure pw ensures the incompressibility with the weak equation:

The isochoric strain energy density is defined by a function of the form

(1)

The three terms on the right-hand side of Equation 1 depend on invariants of the right Cauchy-Green tensor.

The first term describes the mechanical behavior of the elastic ground substance. The isotropic function W1 depends on one material parameter and the first isochoric invariant

, defined in the same fashion as for a neo-Hookean material (see Ref. 3 for more details on this invariant)

(2)

The second and third terms on the right-hand side of Equation 1 describe the mechanical contribution of the collagen fiber network. Following Ref. 1, these expressions are written as

(3)

(4)

here, the fiber network is reduced to two families of fibers with material properties k1 and k2.

The deformation of each fiber family is measured by the invariants I4 and I6. You can find a detailed background for this formulation in Ref. 1 and Ref. 2.

Briefly, a family i of fibers is defined by a vector field a0i in the undeformed direction. The fibers deform under the action of the isochoric deformation gradient, so that

is the deformed fiber configuration. The length of

is the fiber stretch, to be used for the constitutive equations.

The HGO model uses the square of the fiber stretches computed according to the invariants I4 and I6

(5)

(6)

Also, the angle β is the relative angle between a01 and a02.

The mechanical properties of both the media and adventitia are governed by these expressions. Each layer has a distinct set of material parameters c, k1, k2, and the initial fiber directions a01 and a02 are aligned at different angles, as shown in Figure 2.

Figure 2: Angle between the two fiber families in the adventitia.

Material

The material parameters for the media and the adventitia are given in the following table.


Material properties (Ref. 1)	Value, media	Value, adventitia
c	3[kPa]	0.3[kPa]
k1	2.3632[kPa]	0.5620[kPa]
k2	0.8393	0.7112
β	29[deg]	62[deg]

Results and Discussion

The model computes the static response to the applied boundary conditions. Figure 3 displays the fiber layout of both the media and the adventitia fiber families.

Figure 3: Fiber layout in the undeformed configuration of the media (inner, red) and the adventitia (outer, blue). Note the different angles between fiber families.

Figure 4 shows the radial stress distribution through the thickness of the wall at an axial stretch of 1.9 and an internal pressure of 160 mmHg.

Figure 4: Radial stress distribution in the artery wall at an axial stretch of 1.9 and 160 mmHg internal pressure.

Figure 5 plots the internal pressure against the expansion of the inner radius for the entire load range. The results are in good agreements with the data reproduced from Ref. 1.

Figure 5: Plot of internal pressure vs. inner radius for three different axial stretches. Data reproduced from Ref. 1 (circles) coincides well with the model results.

References

1. G. Holzapfel, T. Gasser, and R. Ogden, “A New Constitutive Framework for Arterial Wall Mechanics and a Comparative Study of Material Models,” J. Elasticity, vol. 61, pp. 1–48, 2000.

2. G. Holzapfel, Nonlinear Solid Mechanics: A Continuum Approach for Engineering, John Wiley & Sons, 2000.

3. Nonlinear Structural Materials Module User’s Guide, COMSOL Multiphysics.

Notes About the COMSOL Implementation

The most important aspect in this model is the implementation of the HGO material model with a user-defined strain energy density function. The initial fiber directions are user-defined variables.

There are two fiber families in each arterial layer. The mathematical expressions in the HGO model are the same for the media as for the adventitia, except for the different material parameters. Only the expressions for the media are discussed below, denoted by the index M. Replace this index with the index A to obtain the expressions for the adventitia.

The vector field a01M of the first fiber family in the media is represented by the components a01M1 (radial), a01M2 (azimuthal), and a01M3 (axial). The second fiber family a02M in the media has the components a02M1, a02M2, and a02M3.

Each fiber family has an associated invariant computed according to Equation 5. The internal variables for the isochoric right Cauchy-Green deformation tensor

in the local coordinate system are solid.CIelIJ (Ref. 3), where I and J are indexes running from 1 to 3. In the , 1 represents the radial direction, 2 the azimuthal and 3 the axial direction.

In the settings for the hyperelastic material, select the cylindrical system as a reference coordinate system. The invariant associated with the first fiber family is named I4CIel, and according to Equation 5 it is computed as

a01M1*solid.CIel11*a01M1+2*a01M1*solid.CIel12*a01M2+

2*a01M1*solid.CIel13*a01M3+a01M2*solid.CIel22*a01M2+

2*a01M2*solid.CIel23*a01M3+a01M3*solid.CIel33*a01M3.

The invariant I6CIel is defined similarly, but it uses the components of the second fiber family, a02M1, a02M2, and a02M3.

Once the fiber directions and related invariants are defined, implement the strain energy density functions. Equation 2 defines an isotropic isochoric function. The corresponding invariant

is defined by the variable solid.I1CIel, thus the expression in the media reads cM/2*(solid.I1CIel-3).

The second and third terms on the right-hand side of Equation 1 are the anisotropic strain energy functions for the fiber families a01M and a02M. With the definition of the fiber invariants in the media for the fiber family a01M, Equation 3 becomes

k1M/(2*k2M)*(exp(k2M*(I4CIel-1)^2)-1)*(I4CIel>1).

The last factor, (I4CIel>1), evaluates to zero if the fiber stretch is smaller than one. This means that the fibers only contribute to tensile stress. To implement Equation 4, simply replace I4CIel with I6CIel, which corresponds to the second fiber family a02M.

This model also shows how to use the Curvilinear Coordinates interface to plot the configuration of the fiber families. The interface takes care of the mapping from the global coordinate system in 2D axisymmetric component to the cylindrical coordinates system used in the revolved result plots.

Nonlinear_Structural_Materials_Module/Hyperelasticity/arterial_wall_mechanics

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Use Curvilinear Coordinates to visualize the fiber layout. Add it four times, once for each fiber family.

In the tree, select .

Click four times.

Click

In the tree, select .

Click

Global Definitions

Load all model parameters from a file containing parameters for the geometry, the material properties and the boundary conditions.

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file arterial_wall_mechanics_parameters.txt.

Now add an interpolation function for importing the pressure versus radius data reproduced from Ref. 1. Use it for comparison.

Interpolation 1 (int1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click .

Browse to the model’s Application Libraries folder and double-click the file arterial_wall_mechanics_pressure_radius.txt.

This file contains the data adapted from Ref. 1.

In the text field, type 1.

Click .

Find the subsection. In the table, enter the following settings:


Function name	Position in file
hgo_pr_1_5	1
hgo_pr_1_7	2
hgo_pr_1_9	3

Locate the section. In the text field, type kPa.

In the text field, type mm.

Geometry 1

Construct the model geometry by first drawing two circular sections on a work plane. Then, form a difference between them and extrude it.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type Ro-Ri.

In the text field, type L.

Locate the section. In the text field, type Ri.

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (mm)
Layer 1	HA

Clear the check box.

Select the check box.

In the toolbar, click

Definitions

Load all model variables from files. These define the initial directions of all fiber families. The files also contain the expressions for the user defined strain energy density functions.

Fiber Directions, Media

In the toolbar, click

and choose .

In the window for , type Fiber Directions, Media in the text field.

Locate the section. From the list, choose .

Select Domain 1 only.

Locate the section. Click

Browse to the model’s Application Libraries folder and double-click the file arterial_wall_mechanics_fiber_directions_media.txt.

Fiber Directions, Adventitia

In the toolbar, click

and choose .

In the window for , type Fiber Directions, Adventitia in the text field.

Locate the section. From the list, choose .

Select Domain 2 only.

Locate the section. Click

Browse to the model’s Application Libraries folder and double-click the file arterial_wall_mechanics_fiber_directions_adventitia.txt.

Strain Energy Density, Media

In the toolbar, click

and choose .

In the window for , type Strain Energy Density, Media in the text field.

Locate the section. From the list, choose .

Select Domain 1 only.

Locate the section. Click

Browse to the model’s Application Libraries folder and double-click the file arterial_wall_mechanics_strain_energy_density_media.txt.

Strain Energy Density, Adventitia

In the toolbar, click

and choose .

In the window for , type Strain Energy Density, Adventitia in the text field.

Locate the section. From the list, choose .

Select Domain 2 only.

Locate the section. Click

Browse to the model’s Application Libraries folder and double-click the file arterial_wall_mechanics_strain_energy_density_adventitia.txt.

Solid Mechanics (solid)

Hyperelastic Material 1

In the window, under right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

From the list, choose .

The HGO model is incompressible, select to apply a mixed formulation.

In the Wsiso text field, type W1+W4+W6.

This is the isochoric part of the strain energy density function. The functions W1, W4 and W6 are defined with different properties in the media and adventitia.

From the ρ list, choose . In the associated text field, type 1100.

In the window, click .

In the window for , locate the section.

From the list, choose .

Symmetry Plane 1

In the toolbar, click

and choose .

Select Boundaries 2 and 5 only.

Prescribed Displacement 1

In the toolbar, click

and choose .

Select Boundaries 3 and 6 only.

In the window for , locate the section.

Select the check box.

In the u 0 z text field, type (lambda_z-1)*L.

Use lambda_z as a continuation parameter in the solver to vary axial stretch.

Boundary Load 1

In the toolbar, click

and choose .

Select Boundary 1 only.

In the window for , locate the section.

From the list, choose .

In the p text field, type p_i.

Use p_i as a continuation parameter in the solver settings to vary the internal pressure.

Curvilinear Coordinates (cc)

Now set up the curvilinear coordinates for all fiber families by adding a user defined vector field in the cylindrical system for the components of the fiber directions.

In the window, under click .

Select Domain 1 only.

User Defined 1

In the toolbar, click

and choose .

Set the components of the vector field to the cylindrical components of the first fiber family.

In the window for , locate the section.

Specify the u vector as


a01M1	R
a01M2	PHI
a01M3	Z

Set up the other curvilinear coordinates similarly according to the table:


	Curvilinear Coordinates (cc)	Curvilinear Coordinates 2 (cc2)	Curvilinear Coordinates 3 (cc3)	Curvilinear Coordinates 4 (cc4)
Domain selection	1	1	2	2
u vector, R component	a01M1	a02M1	a01A1	a02A1
u vector, PHI component	a01M2	a02M2	a01A2	a02A2
u vector, Z component	a01M3	a02M3	a01A3	a02A3

Mesh 1

Mapped 1

In the toolbar, click

Distribution 1

Right-click and choose .

Select Boundary 2 only.

In the window for , locate the section.

In the text field, type 6.

Distribution 2

In the window, right-click and choose .

Select Boundary 5 only.

In the window for , locate the section.

In the text field, type 4.

In the window, right-click and choose .

Study 1

Now set up a study to compute the static response of the artery segment subject to combined axial stretch and internal pressure.

In the window, click .

In the window for , locate the section.

Clear the check box.

You will not need the default plots in this model.

Step 1: Stationary

In the window, under click .

In the window for , click to expand the section.

Select the check box.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
lambda_z (Axial stretch)	1.5 1.7 1.9

The parameter lambda_z controls the axial stretch.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
p_i (Internal pressure)	range(0,5,160)	mmHg

Use p_i to vary the internal pressure from 0 to 160 mmHg with steps of 5 mmHg.

From the list, choose .

Using the option for is suitable for this kind of multiparameter sweep with continuation.

Solution 1 (sol1)

In the toolbar, click

Using constant prediction for the continuation sweep improves convergence when the solution is very nonlinear in the swept parameter.

In the window, expand the node.

In the window, expand the node, then click .

In the window for , click to expand the section.

From the list, choose .

In the toolbar, click

Results

In the window, expand the node.

Before you examine the results, create datasets. Those datasets create the cylindrical geometry from the 2D axisymmetric plane that was used for the computation.

Study 1/Solution 1 (sol1)