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 tissue 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 deformations.

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

This example 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 a 2D axisymmetric model. The 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 following set of boundary conditions replicate these experiments: On the bottom surface, a symmetry boundary condition is applied, which allows the bottom end of the artery to expand freely in the radial direction. On the top surface, prescribed displacements in the axial direction account for the axial stretching. The 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 described mathematically 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 elastic, isochoric right Cauchy–Green deformation tensor. Here, the first term describes the mechanical behavior of the isotropic ground substance. The strain energy density function W1 depends on one material parameter c and the first invariant

, defined in the same fashion as for an incompressible neo-Hookean material (see Ref. 3 for more details)

(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. 4, these expressions are written as

(3)

(4)

Herein, the fiber network is reduced to two families of fibers with equal material properties k1, k2, and k3. The parameter k1 represents the fiber stiffness (SI unit: Pa), k2 is a dimensionless tuning parameter, and k3 is the fiber dispersion. A family i of fibers is defined by a vector field a0i of unit length in the reference configuration. The fibers deform under the action of the isochoric deformation gradient, so that

describes the fiber vector in the current configuration. The length of

is the fiber stretch, which can be used in the constitutive equations.The HGO model uses the square of the fiber stretches according to the invariants

(5)

Following the example in Ref. 1, the two fiber families are assumed to be oriented at an angle

with respect to the circumferential direction of the artery. The subscript

indicates that the angle differs between the media and the adventitia. You can find a detailed background of the HGO model formulation in Ref. 1, Ref. 2, and Ref. 4.

In this example, the mechanical properties of both the media and the adventitia are assumed to be governed by these expressions. Each layer has a distinct set of material parameters c, k1, and k2, and k3, and the initial fiber directions a0i are aligned at different angles.

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 initial alignment of the fiber families in the media and the adventitia is shown in Figure 2.

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

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

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

The internal pressure is plotted as a function of the inner radius of the artery for pressures from 0 to 160 mmHg and axial stretches of 1.5, 1.7, and 1.9 in Figure 4. The results are in excellent agreement with the data reproduced from Ref. 1.

Figure 4: Plot of internal pressure vs. inner radius for three different axial stretches. Data reproduced from Ref. 1 (circles) coincide 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.

4. T. Gasser, R. Ogden, and G. Holzapfel, “Hyperelastic modelling of arterial layers with distributed collagen fibre orientations”, J. R. Soc. Interface, vol 3, pp. 15–35, 2006.

Notes About the COMSOL Implementation

The most important aspect in this model is the implementation of the HGO material model with the feature.

There are two fiber families in each arterial layer. The mathematical expressions in the HGO model are the same for the media (M) and 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 is aligned with the fiber directions.

The mechanical behavior of the elastic ground substance is defined in the parent 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 feature is used to compute such strain energy densities. Add a feature for each fiber family, 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 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 second axis (x2).

The feature automatically computes the corresponding invariant

according to Equation 5 in order to compute the strain energy function Wi in Equation 3.

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

Activate the option Contribute to total stress to add the fiber 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 orientation of the fiber families (see Figure 2). The COMSOL implementation takes care of the mapping from the global coordinate system in a 2D axisymmetric model to the cylindrical coordinate 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 drawing a circumferential cross-section of the artery. Use to divide it into two domains corresponding to the Media and the Adventitia.

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 model parameters loaded.

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 in 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 the convergence when the solution is very nonlinear in the sweep 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. These datasets create the 3D 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 .

In the window for , type Sector Revolution in the text field.

Click to expand the section. In the text field, type -90.

In the text field, type 225.

Full Revolution

In the toolbar, click

and choose .

In the window for , type Full Revolution in the text field.

Now duplicate the datasets and add a selection for the media. Use these in one of the plots below.

Media

In the window, under right-click and choose .

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

Selection

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Domain 1 only.

Media Sector Revolution

In the toolbar, click

and choose .

In the window for , type Media Sector Revolution in the text field.

Locate the section. From the list, choose .

Locate the section. In the text field, type -40.

In the text field, type 140.

Click to expand the section.

Now duplicate the datasets and add a selection for the adventitia. Use these in one of the plots below.

Adventitia

Right-click and choose .

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

Selection

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Domain 2 only.

Adventitia Sector Revolution

In the toolbar, click

and choose .

In the window for , type Adventitia Sector Revolution in the text field.

Locate the section. From the list, choose .

Locate the section. In the text field, type -50.

In the text field, type 160.

Create a 3D plot group for the radial stress distribution.

Radial Stress

In the toolbar, click

In the window for , type Radial Stress in the text field.

Surface 1

In the toolbar, click

In the window for , click in the upper-right corner of the section. From the menu, choose .

Locate the section. From the list, choose .

Deformation 1

In the toolbar, click

In the window for , locate the section.

Select the check box. In the associated text field, type 1.

In the toolbar, click

Click the

button in the toolbar.

Create a 1D plot group to compare the pressure versus radius relationship with the data from Ref. 1.

Pressure vs. Radius

In the toolbar, click

and choose .

In the window for , type Pressure vs. Radius in the text field.

Click to expand the section. From the list, choose .

Locate the section.

Select the check box. In the associated text field, type Inner radius (mm).

Select the check box. In the associated text field, type Internal pressure (mmHg).

Locate the section. From the list, choose .

Point Graph 1

Right-click and choose .

Select Point 1 only.

In the window for , locate the section.

In the text field, type p_i.

In the field, type mmHg.

Click in the upper-right corner of the section. From the menu, choose .

In the toolbar, click

Click to expand the section. Select the check box.

Find the subsection. Clear the check box.

Global 1

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the table, enter the following settings:


Expression	Unit	Description
p_i	mmHg	Internal pressure

Locate the section. From the list, choose .

In the text field, type hgo_pr_1_5(p_i).

This is the interpolation function with data from Ref. 1. It returns the inner radius as a function of the internal pressure at an axial stretch of 1.5.