PDF
Arterial Wall Viscoelasticity
Introduction
Arteries are blood vessels that carry freshly oxygenated blood from the heart throughout the rest of the body. The Holzapfel–Gasser–Ogden (HGO) model proposes a mechanical description of the young healthy arteries in Ref. 1 based on anisotropic hyperelastic properties, which is implemented in the Arterial Wall Mechanics example. Here we study the dynamic behavior of the artery, based on Ref. 2 and calculate the time-dependent response given a sudden axial stretching.
Model Definition
The geometry, physics interface, and material models are the same as in the example Arterial Wall Mechanics, see also Figure 1.
Figure 1: Carotid artery section made of a media layer and an adventitia layer.
The two modeled layers (media and adventitia) are described by the anisotropic HGO hyperelastic material model:
Only the media layer includes viscoelastic behavior. The generalized Maxwell viscoelastic model is used to represent relaxations time at different time-scales (Ref. 2):
(1)
For each branch of the generalized Maxwell model, the viscoelastic stress follows:
where τm is the relaxation time and βm is the energy factor per branch. Applying a variable change qm = βmSiso − Qm gives
(2)
In this example, a generalized Maxwell viscoelastic model with five branches is used with the following values taken from Ref. 2:
Branch
Energy factor
Relaxation time
1
0.3353
0.001 s
2
0.286
0.01 s
3
0.298
0.1 s
4
0.285
1 s
5
0.348
10 s
The artery is first loaded with an internal pressure of 100 mmHg and an initial axial stretch of 1.5. After initialization, the stretch is increased to 1.7 and the viscoelastic relaxation is calculated.
Results and Discussion
The total force is computed by integrating the axial stress on the top section surface. The plot shown in Figure 2 is similar to the force relaxation presented in Ref. 2. The force relaxes almost linearly from 10-3 s to 10 s due to the wide range of relaxation times.
Moreover, Figure 3 shows the evolution of total viscoelastic stress and viscoelastic stress for each branch over time. The viscoelastic stress for the branches is calculated with the expression Qm = βmSiso − qm. The viscoelastic stress relaxes from its initial value to zero.
Figure 2: Relaxation of the axial force after stretching.
Figure 3: Variation of viscoelastic stress. Both the total stress and the stress in each branch is shown.
Notes About the COMSOL Implementation
A stationary study step is needed to prestress the artery with initial pressure and stretch. As this initial state is assumed to be a steady-state, the static stiffness in Viscoelasticity node must be set to Long-term. This ensures that the viscoelastic model has no effects in the stationary step.
References
1. G.A. Holzapfel, T.C. Gasser, and R.W. 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.A. Holzapfel, T.C. Gasser, and M. Stadler, “A structural model for the viscoelastic behavior of arterial walls: Continuum formulation and finite element analysis,” Eur. J. Mech. A Solids, vol. 21, no. 3, pp. 441–463, 2002
Application Library path: Nonlinear_Structural_Materials_Module/Viscoelasticity/arterial_wall_viscoelasticity
Modeling Instructions
Application Libraries
1
From the File menu, choose Application Libraries.
2
In the Application Libraries window, select Nonlinear Structural Materials Module > Hyperelasticity > arterial_wall_mechanics in the tree.
3
Click  Open.
Global Definitions
Parameters 1
Set the parameters used for loads and constraints in the time dependent study. Parameter t is the time, and is needed to define the prescribed displacement in the stationary study step.
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
In the table, enter the following settings:
Name
Expression
Value
Description
p_i
100[mmHg]
13332 Pa
Internal pressure
lambda_z
1.7
1.7
Axial stretch
lambda_z0
1.5
1.5
Initial axial stretch
t
0[s]
0 s
Time
Create a step function to apply the stretch in the time-dependent study.
Step 1 (step1)
1
In the Home toolbar, click  Functions and choose Global > Step.
2
In the Settings window for Step, type lambda in the Function name text field.
3
Locate the Parameters section. In the Location text field, type 5e-6[s].
4
In the From text field, type lambda_z0.
5
In the To text field, type lambda_z.
6
Click to expand the Smoothing section. In the Size of transition zone text field, type 1e-5.
7
Click  Plot.
Component 1 (comp1)
In the Model Builder window, expand the Component 1 (comp1) node.
Solid Mechanics (solid)
Hyperelastic Material (Media)
Add a generalized Maxwell viscoelasticity model according to Ref. 1.
1
In the Model Builder window, expand the Component 1 (comp1) > Solid Mechanics (solid) node, then click Hyperelastic Material (Media).
Viscoelasticity 1
1
In the Physics toolbar, click  Attributes and choose Viscoelasticity.
2
Select Domain 1 only.
3
In the Settings window for Viscoelasticity, locate the Viscoelasticity Model section.
4
Click  Add four times.
5
In the table, enter the following settings:
Branch
Energy factor (1)
Relaxation time (s)
1
0.353
0.001
2
0.286
0.01
3
0.298
0.1
4
0.285
1
5
0.348
10
Prescribed Displacement 2
1
In the Model Builder window, under Component 1 (comp1) > Solid Mechanics (solid) right-click Prescribed Displacement 1 and choose Duplicate.
2
In the Settings window for Prescribed Displacement, locate the Prescribed Displacement section.
3
In the u 0 z text field, type (lambda(t)-1)*L.
Definitions
Integration 1 (intop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
In the Settings window for Integration, locate the Source Selection section.
3
From the Geometric entity level list, choose Boundary.
4
Select Boundaries 3 and 6 only.
Mesh 2
In the Mesh toolbar, click Add Mesh and choose Add Mesh.
Mapped 1
In the Mesh toolbar, click  Mapped.
Distribution 1
1
In the Model Builder window, right-click Mapped 1 and choose Distribution.
2
Select Boundaries 2 and 3 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
From the Distribution type list, choose Predefined.
5
In the Number of elements text field, type 20.
6
In the Element ratio text field, type 10.
7
Select the Reverse direction checkbox.
Distribution 2
1
In the Model Builder window, right-click Mapped 1 and choose Distribution.
2
Select Boundary 5 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
In the Number of elements text field, type 8.
5
In the Model Builder window, right-click Mesh 2 and choose Build All.
Add a new study. The first stationary step is used to prestress the artery with a stretch and an internal pressure. Select Long-term stiffness in the Viscoelasticity node, so the viscoelastic effect is disabled. The time-dependent step is used to compute the dynamic response to an additional stretch.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies > Stationary.
4
Click the Add Study button in the window toolbar.
5
In the Home toolbar, click  Add Study to close the Add Study window.
Study 2
Step 1: Stationary
1
In the Settings window for Stationary, click to expand the Study Extensions section.
2
Select the Auxiliary sweep checkbox.
3
Click  Add.
4
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
p_i (Internal pressure)
range(0,10,100)
mmHg
Step 2: Time Dependent
1
In the Study toolbar, click  Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type range(0,0.5e-6,9.5e-6) 10^{range(-5,0.2,2.4)}.
Solution 2 (sol2)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 2 (sol2) node, then click Time-Dependent Solver 1.
3
In the Settings window for Time-Dependent Solver, click to expand the Time Stepping section.
4
From the Steps taken by solver list, choose Intermediate.
5
Find the Algebraic variable settings subsection. From the Error estimation list, choose Exclude algebraic.
6
In the Model Builder window, click Study 2.
7
In the Settings window for Study, locate the Study Settings section.
8
Clear the Generate default plots checkbox.
9
In the Study toolbar, click  Compute.
Results
Integrate the axial stress on the top surfaces to calculate the reaction force and reproduce Figure 2.
Force Relaxation
1
In the Model Builder window, expand the Results > Datasets node.
2
Right-click Results and choose 1D Plot Group.
3
In the Settings window for 1D Plot Group, type Force Relaxation in the Label text field.
4
Locate the Data section. From the Dataset list, choose Study 2/Solution 2 (sol2).
5
From the Time selection list, choose Interpolated.
6
In the Times (s) text field, type 10^{range(-5,0.2,2.4)}.
7
Click the  x-Axis Log Scale button in the Graphics toolbar.
8
In the Force Relaxation toolbar, click  Plot.
Global 1
1
Right-click Force Relaxation and choose Global.
2
In the Settings window for Global, locate the y-Axis Data section.
3
In the table, enter the following settings:
Expression
Unit
Description
intop1(solid.sz)
mN
Axial force
4
In the Force Relaxation toolbar, click  Plot.
Add a new plot group to plot viscoelastic stress and reproduce Figure 3.
Viscoelastic Stress
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Viscoelastic Stress in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 2/Solution 2 (sol2).
4
From the Time selection list, choose Interpolated.
5
In the Times (s) text field, type 10^{range(-5,0.2,2.4)}.
6
Click to expand the Title section. From the Title type list, choose Manual.
7
In the Title text area, type Radial viscoelastic stress.
8
Locate the Plot Settings section.
9
Select the y-axis label checkbox. In the associated text field, type Viscoelastic Stress (N/m<sup>2</sup>).
10
Locate the Legend section. From the Position list, choose Lower right.
Point Graph 1
1
Right-click Viscoelastic Stress and choose Point Graph.
2
Select Point 1 only.
3
In the Settings window for Point Graph, click Replace Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1) > Solid Mechanics > Stress > Viscoelastic stress tensor (material and geometry frames) - N/m² > solid.SqGpRR - Viscoelastic stress tensor, RR-component.
4
Click to expand the Legends section. Select the Show legends checkbox.
5
From the Legends list, choose Manual.
6
In the table, enter the following settings:
Legends
Total viscoelastic stress
Point Graph 2
1
Right-click Point Graph 1 and choose Duplicate.
2
In the Settings window for Point Graph, locate the y-Axis Data section.
3
In the Expression text field, type solid.gpeval(solid.hmm1.vis1.betavm1*solid.Sliso11-solid.hmm1.vis1.qm1_11).
4
Locate the Legends section. In the table, enter the following settings:
Legends
Branch 1
5
Duplicate the point graph four times and change properties according to the table:
Point graph
Expression of y-Axis Data
Legends
Point Graph 3
solid.gpeval(solid.hmm1.vis1.betavm2*solid.Sliso11-solid.hmm1.vis1.qm2_11)
Branch 2
Point Graph 4
solid.gpeval(solid.hmm1.vis1.betavm3*solid.Sliso11-solid.hmm1.vis1.qm3_11)
Branch 3
Point Graph 5
solid.gpeval(solid.hmm1.vis1.betavm4*solid.Sliso11-solid.hmm1.vis1.qm4_11)
Branch 4
Point Graph 6
solid.gpeval(solid.hmm1.vis1.betavm5*solid.Sliso11-solid.hmm1.vis1.qm5_11)
Branch 5
Point Graph 1
1
In the Model Builder window, click Point Graph 1.
2
In the Settings window for Point Graph, click to expand the Coloring and Style section.
3
Find the Line markers subsection. From the Marker list, choose Plus sign.
4
From the Positioning list, choose Interpolated.
Viscoelastic Stress
1
Click the  x-Axis Log Scale button in the Graphics toolbar.
2
In the Model Builder window, click Viscoelastic Stress.
3
In the Viscoelastic Stress toolbar, click  Plot.
You can group plot groups to improve the clarity of the Results tree.
Radial Stress
In the Model Builder window, under Results right-click Radial Stress and choose Group.
Stationary Results
In the Settings window for Group, type Stationary Results in the Label text field.
Pressure vs. Radius
In the Model Builder window, right-click Pressure vs. Radius and choose Stationary Results.
Stationary Results
In the Model Builder window, right-click Stationary Results and choose Move Down.
Force Relaxation
In the Model Builder window, under Results right-click Force Relaxation and choose Group.
Viscoelasticity Results
In the Settings window for Group, type Viscoelasticity Results in the Label text field.
Viscoelastic Stress
In the Model Builder window, right-click Viscoelastic Stress and choose Viscoelasticity Results.
Disable the new prescribed displacement and the viscoelasticity nodes in study 1 in order to keep it in its original state.
Study 1
Step 1: Stationary
1
In the Model Builder window, expand the Study 1 node, then click Step 1: Stationary.
2
In the Settings window for Stationary, locate the Physics and Variables Selection section.
3
Select the Modify model configuration for study step checkbox.
4
In the tree, select Component 1 (comp1) > Solid Mechanics (solid), Controls spatial frame > Prescribed Displacement 2.
5
Click  Disable.
6
In the tree, select Component 1 (comp1) > Solid Mechanics (solid), Controls spatial frame > Hyperelastic Material (Media) > Viscoelasticity 1.
7
Click  Disable.