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 the Fiber feature.

There are two fiber families in each arterial layer. The mathematical expressions in the HGO model are the same for the media (M) as for the adventitia (A), except for the different material parameters.

The initial fiber directions are identified with the use of user defined rotated coordinate systems. Four coordinate systems are necessary, one for each fiber family. They are oriented in such a way that their second axis are aligned with the fiber directions (a01M, a02M, a01A, a02A).

The mechanical behavior of the elastic ground substance is defined in the Hypereleastic node (one for the media and one for the adventitia) with the use of an incompressible Neo-Hookean model to define the isotropic function W1 as in Equation 2.

The mechanical contribution of the collagen fiber network is accounted for with an anisotropic contribution to the isochoric strain energy (W4+W6). The Fiber features are used to computed such strain energy densities. Add a Fiber feature for each family of fiber, two for the media and two for the adventitia. Only the operation for the first family of fibers in the media is discussed below. Similar operations are performed for all other fiber families.

In the settings for the Fiber feature, select the appropriate rotated reference system as a reference coordinate system. Then in the Orientation section select the orientation of the fiber (a01M) to be aligned with the x2-axis.

The Fiber feature automatically compute the corresponding invariants (I4M) as in Equation 5 in order to compute the strain energy function W4 as in Equation 3.

The option Stiffness in tension only is already activated for the HGO model and it evaluates the fibers strain energy to zero if the fiber stretch is smaller than one. This means that the fibers only contribute to tensile stress.

Activate the option Contribute to total stress to add the fibers strain energy directly to the total energy. This option allows to consider the fibers and the ground substance as a single anisotropic material and the resulting stress will account for all contributions.

This model also shows how to use the fiber direction variables 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 .

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 table, enter the following settings:


Argument	Unit
Column 1	kPa

In the table, enter the following settings:


Function	Unit
hgo_pr_1_5	mm
hgo_pr_1_7	mm
hgo_pr_1_9	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

Define coordinate systems. These are used to define the initial directions of all fiber families.

Rotated System Media Fiber Family 1

In the toolbar, click

and choose .

In the window for , type Rotated System Media Fiber Family 1 in the text field.

Locate the section. From the list, choose .

Find the subsection. In the β text field, type betaM.

Locate the section. From the list, choose .

Rotated System Media Fiber Family 2

Right-click and choose .

In the window for , type Rotated System Media Fiber Family 2 in the text field.

Locate the section. Find the subsection. In the β text field, type -betaM.

Rotated System Adventitia Fiber Family 1

Right-click and choose .

In the window for , type Rotated System Adventitia Fiber Family 1 in the text field.

Locate the section. Find the subsection. In the β text field, type betaA.

Rotated System Adventitia Fiber Family 2

Right-click and choose .

In the window for , type Rotated System Adventitia Fiber Family 2 in the text field.

Locate the section. Find the subsection. In the β text field, type -betaA.

Materials

Create the materials for the media and the adventitia.

Material (Media)

In the window, under right-click and choose .

In the window for , locate the section.

In the list, select .

Click

Select Domain 1 only.

In the text field, type Material (Media).

Material (Adventitia)

Right-click and choose .

In the window for , type Material (Adventitia) in the text field.

Select Domain 2 only.

Solid Mechanics (solid)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Hyperelastic Material (Media)

In the toolbar, click

and choose .

In the window for , type Hyperelastic Material (Media) in the text field.

Select Domain 1 only.

Hyperelastic Material (Adventita)

In the toolbar, click

and choose .

In the window for , type Hyperelastic Material (Adventita) in the text field.

Select Domain 2 only.

Locate the section. From the list, choose .

Hyperelastic Material (Media)

In the window, click .

In the window for , locate the section.

From the list, choose .

Fiber Family 1

In the toolbar, click

and choose .

In the window for , type Fiber Family 1 in the text field.

Locate the section. From the list, choose .

Locate the section. From the a list, choose .

Locate the section. Select the check box.

Fiber Family 2

Right-click and choose .

In the window for , type Fiber Family 2 in the text field.

Locate the section. From the list, choose .

Hyperelastic Material (Adventita)

In the window, under click .

Fiber Family 1

In the toolbar, click

and choose .

In the window for , type Fiber Family 1 in the text field.

Locate the section. From the list, choose .

Locate the section. Select the check box.

Locate the section. From the a list, choose .

Fiber Family 2

Right-click and choose .

In the window for , type Fiber Family 2 in the text field.

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.

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.

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 .

Materials

Material (Media) (mat1)

Before solving all necessary materials properties should be defined using the loaded model parameters.

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Lamé parameter μ	muLame	cM	N/m²	Lamé parameters
Density	rho	rhoM	kg/m³	Basic
Fiber stiffness	k1HGO	k1M	Pa	Holzapfel-Gasser-Ogden
Model parameter	k2HGO	k2M	1	Holzapfel-Gasser-Ogden
Fiber dispersion	k3HGO	k3M	1	Holzapfel-Gasser-Ogden

Material (Adventitia) (mat2)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Lamé parameter μ	muLame	cA	N/m²	Lamé parameters
Density	rho	rhoA	kg/m³	Basic
Fiber stiffness	k1HGO	k1A	Pa	Holzapfel-Gasser-Ogden
Model parameter	k2HGO	k2A	1	Holzapfel-Gasser-Ogden
Fiber dispersion	k3HGO	k3A	1	Holzapfel-Gasser-Ogden

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.

Sector Revolution

In the window, expand the node.

Right-click and choose .