Layered Vascular Vessel with Yeoh Constitutive Model

IntroductionComparison
Analytical solution approximationNomenclature
Triple layered modelReferences
Homogenized models

Introduction

This models analyses the stress state of a coronary vessel under physiological conditions. The model makes use of a Yeoh hyperelastic material model with model parameters extracted from experimental procedures.

The coronary arteries are major blood vessels supplying blood to the heart. There are several diseases linked to the coronary arteries that can compromise an entire organism. There are several diseases related to coronary arteries. A common type of heart disease is caused by plaque buildup in the wall of coronary arteries. Also, a coronary aneurism represent a dangerous pathology. It is defined as a dilatation of the coronary arteries in which the wall exceeds the normal diameter. These diseases are strongly linked to the wall of the arteries. The artery walls are an active structure that continuously remodels during the entire life of a human. It senses the bio-chemo-mechanical environment to guarantee homeostasis and the best blood supply to heart tissues [Mishani, 2021].

1.gif

This figure shows the complexity reduction made from the anatomical structure of the coronary vessel to a simplified model that maintains mechanical properties of interest. The 2D structure is made up of three layers and symmetry constraints (only sliding in the edge-parallel direction) which allows representing the entire annulus.

The structure of the wall of an artery is complex. From a topological point of view, it consists of three different layers. The inner one, called tunica intima, is a layer made of epithelium and connective tissues with elastic fibers. The middle layer, the tunica media, is primarily made of smooth muscle cells and can regulate the arterial diameter due to local and central stimuli. The outer layer, the tunica adventitia, is made of connective tissue with collagen and elastic fibers. It links the artery to the surrounding structures.

Here a comparison between a triple-layered vessel and a homogenous one is presented. The three layers represent the tunica intima, media and adventitia. Each one with different mechanical properties from the literature. Geometry is reduced to three concentric cylinders with a diameter representative of the anatomy of the main left coronary artery [Gholipour, 2018]. Furthermore, two different techniques to homogenize the vessel properties towards a homogeneous single layer structure are implemented.

Four different modeling approaches are shown: first a simple analytical model is presented. Then, in a second model, the usage of the three layer model with different material parameters in each layer to represent the tunica intima, media and adventitia is implemented. Then, the next two models represent simplified homogenisations. One of the homogenisation approaches is to use the material parameters of the inner layer, the tunica intima layer and in the second homogenisation approach an area average of all three material layers is used.

In order to show the difference between a multi layer constitutive model and a single layer model an analysis of a uniaxial tensile traction test is performed. This is done by applying a stretch along a single direction and it leads to and λ2=λ3.

To enforce incompressibility λ1λ2λ3=1 is set. It is then possible to express the first strain invariant I1 with respect to the uniaxial stretch, as:

And the Cauchy stress over the first principal direction reads

with the pressure

such that

Specify different material parameters for the different layers:

Then it is possible to plot the material constitutive models as bi-dimensional chart, strain vs stress. With biological tissue, it is often useful work with kPa and a scale factor is introduced to rescale quantities.

Introduce a scale factor useful to rescale quantities:
Define a helper function to compute the uniaxial stress in kPa:
Visualize the uniaxial stress for the three layers:

Analytical solution approximation

This mechanical problem can be also reduced to a characteristic engineering problem where a cylinder or a sphere is subjected to an internal pressure p. This is called a pressured vessel and applications for it arise in many areas.

7.gif

This figure shows an ideal tube with length l, thickness t and internal pressure . It is possible to read circumferential stress and radial stress . Because the of the the open ends of the tube the axial stress can be excluded.

Here the analytical solution is shown for an inflated thicken-wall tube as an upper bound. This can be used to compare the stress state with a FEM solution. With biological tissues, the small strain hypothesis often fails and this leads to several differences between this reduced analytical solution and the numerical one.

The first and the second principal stress directions are shown. The first principal stress is directed across the circumferential direction. The second principal stress direction is along the radial direction.

In the inflated tube, due to the internal pressure, the radial stress is related to this boundary pressure, at least on the internal side. It should also tend to zero on the external side due to the absence of external force or pressure.

Furthermore, this cylinder shows the ratio and this leads to the thick-walled vessels that could be calculated using the Lamè equation. Then it is possible to express the radial and tangential stress as:

where and as

and

with the internal pressure and the external pressure .

The axial stress is close to zero due the absence of axial force or constraints.

Analytical solution of a thick-wall inflated cylinder:
Introduce geometrical quantities:
Visualize the analytical solution:

As it is possible to observe from the plot, the radial stress is related to the boundary conditions: Indeed the radial stress varies from the internal pressure p=140 mmHg 18.6 kPa , on the lumen side of the vessel, to p=0 on the external side. The circumferential stress is strongly related to the properties of the vascular tissue and due to the high non-linearity high differences are present respect to the following FEM analysis.

Triple layered model

The following FEM analysis represent the triple layered model. It uses appropriate simplification in geometry and constitute models in order to reduce computational cost but still obtain accurate results. This relatively simplified analysis could be a starting point to investigate the mechanical environment of arteries and introduce further investigation, eventually with patient-specific geometries, to assess arteries disease risk, development and progression [Holzapfel, 2018].

Load the finite element package and set the $HistoryLength to 0:
Setup model variables:

Each layer is represented as a homogenous material described with the Yeoh hyperelastic model that can reproduce the non-linearity present in the mechanical response of soft tissues such as biological tissue. Each layer uses different material parameters from literature [Gervaso, 2018].

The Yeoh hyperelastic constitutive model is a phenomenological model for representing the deformation of an isotropic, nearly incompressible, non-linear elastic material such as rubber or biological tissues. The model is described by a strain energy density function represented as a power series of the strain invariant.

The Yeoh model implemented is a nearly incompressible model and can be used by specifying the "MaterialMoldel"->"Yeoh" in the parameter pars. The Yeoh model take three parameters called "c1", "c2" and "c3" and it is described from a polynomial function representing the volumetric part of the strain energy density function :

and the first invariant.

For each layer material parameters need to be specified. To make this process easy, a helper function is created to map the properties to each region.

Define a function to map material properties to each layer:
Define marker to specify the different layers:
Specify the model material parameters:

Also, see this note about how to set up of computationally efficient PDE coefficients.

Mesh

In order to reduce the computational cost and better investigate the relationship between each layer, the analysis is reduced to a 2D section of the vessel. Furthermore, due to the rotational symmetry, only a quarter of the annulus is investigated, with corresponding symmetry conditions on the cutting edges.

Generate the boundary for the geometry:
Define marker to specify the different regions:
Generate a mesh:
Visualize the mesh:

It is important to guarantee a sufficient number of elements across the thickness of the layers. At least three layers of elements are needed for each region in order to estimate strain and stress without under-sampling the derivative. More elements would increase accuracy, however the number of elements needs to be balanced with computational cost.

Inspect the number of element for each sub region:

Model setup

For this model, two types of boundary conditions are needed. The symmetry conditions to model the reduced simulation domain and boundary conditions to model internal pressure. For more information on the symmetry condition the section on Displacement Constraints in the Solid Mechanics monograph. With regards to modeling the pressure the goal is to set up a parametric solver with a parametric variable that is the imposed external load. The parameter allows to slowly increase the applied forces, which is needed since the problem is too strongly nonlinear to solve in a single step. For more information on the solving of strongly nonlinear hyperelastic material models see the Neo-Hookean model section in the Hyperelastic Solid Mechanics monograph.

The internal pressure need to be specified as related to the physiological condition of the coronary arteries. A pressure of 140 mmHg [Ramanathan, 2005] is specified.

Specify the total pressure, the maximal value of :

The next step is to specify the positions where this pressure is applied.

Placement of pressure boundary conditions:

The symmetry conditions apply at and , respectively and the displacement is set to 0 in the respective direction which the other direction left free to move. To avoid inaccuracies at the boundary a tolerance is set.

Set a small boundary tolerance:
Create the PDE model with boundary conditions:
Create the parametric solver:

Solution

Due to the high non-linearity of the mechanical problem, the external load is slowly incremented.

Solve for the displacement by slowly increasing the load :

When simulations run for a long time it may make sense to store the obtained solution such that they are readily available later with out re-running the simulation.

Un-comment the following to Export the solution:
Un-comment the following to Import a previously exported solution:

Post processing

After obtaining the solution the data can be processed to investigate stress and strain fields in the domain of the sample.

Stress and strain can be extracted and plotted over both the original and the deformed mesh.

Create and visualize a deformed mesh with the boundary of the original region:

In order to analyze the stress and strain fields, they need to be extracted from the solution.

Compute the strain field from the solution:
Plot the strain field:

The strain fields along the - and -directions are reflected across the bisector at 45°. This is related to the rotational symmetry that this annuls implicitly shows. For this reason, another options is to plot the equivalent strain. The equivalent strain combines all strain components into a scalar and as such is easier to grasp. The equivalent strain is used in the following comparison between homogenous and inhomogenous domains.

Compute the equivalent strain field from the solution:
Plot the equivalent strain field:

It is also interesting to investigate stress distribution across the thickness. The symmetric plots are shown also for the stress across the - and -directions. The Cauchy stress can be plotted showing that it is higher on the inner layer than on the outer.

Compute the Cauchy stress field from the solution:

In order to represent the stress field it is useful to use kPa units, so multiply by a conversion factor.

Plot the Cauchy stress in kPa:

An alternative is to plot the von Mises stress that represents an average across the different components of the stress. Indeed, its plot shows how the stress is higher on the inner layer side and slowly decreases from the inner layer to the outer one.

Compute the von Mises stress:
Compute the von Mises stress on the deformed mesh:
Compute the min and max of the von Mises stress:
Plot the von Mises stress across the deformed vessel:

Also, an investigation of the principal stress direction is useful to understand the stress distribution.

Compute the first and second principal stress directions:
Visualize the first and second principal stress directions:

The principal stresses are eigenvalues of the stress matrix and correspond to the diagonal stress values where the shear stress is 0. For more information on the principal stresses see the section on Principal Directions in the Solid Mechanics monograph. By plotting the principal stress direction it is possible to show how the first one is directed tangentially to the annulus and is related to the stretch across the concentric circumferences. Moreover, the second principal direction appears across the radial direction and the second principal stress is mainly related to the internal boundary pressure and decreases to zero going outward. The following plots show the principal direction while the stress values are shown in the following comparison or explained by the previous analytical solution.

Homogenized models

The previous example made use of a three layer model to describe the artery's wall. The question now is, how does a single layer model compares to the three layer model. In the following, two different homogenization approaches are shown. The first is using the material proprieties of the inner layer, the tunica intima, for all of the simulation domain. This reasonable, especially in pathologies where an intima thickening is present, since in that case the tunica intima is the main load-carrying layer. The second approach uses a homogenization computed as an average of all three layers.

Homogenization with intima values

In order to homogenize the material properties it is possible to use only the properties of the inner layer, the tunica intima. This is also reasonable following previous results showing higher stress in the intima that results in the main load-carrying layer.

Parameters which use the coefficients of the tunica intima:
Create the parametric solver for the intima homogenization:
Solve with the intima value homogenization:
Compute post-processing quantities for the tunica intima homogenized region:

Homogenization with average values

The second homogenization approach uses a mathematical average to homogenize the properties across the three layers into a single one.

Helper function for averaging of material coefficients:

The usage of Round in the above helper function is to avoid numerical instabilities when converting the second Piola stress to a Cauchy stress. Since the average is only an approximation the number is converted to it's nearest integer. Since the integer itself is an exact number the conversion to the Cauchy stress is less likely to run into trouble.

Compute the average of the coefficient and assign it to the new material parameters:
Create the parametric solver for the average homogenization:
Solving with the averaged value homogenization:
Compute post-processing quantities for the averaged value homogenized region:
Compute minimal and maximal strain measures:

Comparison

In this section the three different FEM models are compared.

The layered structure shows more accurate results across the layers interface. However, it requires more than 20 times the computational time of the homogenized ones. Depending on the interest of the simulation it might be reasonable to use one model rather than the other, this requires evaluating the differences in stress and strain results to choose the optimal model.

Compute minimal and maximal equivalent strain between each models:
Plot equivalent strain for the different models:
Compute minimal and maximal von Mises stress:
Compute global maximum and minimum of the von Mises stress to scale color map:
Plot the von Mises stress for the different models:

The homogeneous and homogenized cases are closer together then the layered approach. The difference between the layered and the homogenised models are caused by the presence of an inner material boundary layer for the layered model while no such layer exists in the homogenised models.

Then the stress measures is extracted across the thickness and plot into an uniaxial chart. Due to the symmetry these quantities remain the same for every radial direction.

Create an auxiliary function to highlight different regions in the graphs below:
Create a helper function to extract data on :
Definition of the number of points on which extract the results:
Plot the von Mises stresses across the thickness:
Extract first and second principal stress for the triple-layered model:
Extract first and second principal stress for the homogenous model with the coefficients of the intima:
Extract first and second principal stress for the homogenous model with averaged coefficients:
Plot the principal stresses across the thickness over a radial coordinate:

The first principal stress is related to the stretch in the tangential direction. Indeed, the different models show different strain distributions. These differences are what cause differences in the von Mises stress across the different models as can be seen from the figure.

Plot the principal stresses across the thickness over a radial coordinate:

The second principal stress is similar between each model, indeed this is the stress in the radial direction and in the internal side has to contrast the internal pressure of the boundary value problem (equilibrium law). Then it should go to zero across the thickness due to the absence of external pressure.

It is also interesting that across the interface between the intima and the media layer the stress distribution coincides between the three models. Furthermore, the triple-layered model shows a major load-carrying role for the intima layer with lower stress in the other layers. In the homogeneous and homogenized cases, the stress field continuously decreases from the internal side to the external one.

Nomenclature

References

1.  A. Gholipour, M. H. Ghayesh, A. Zander, R. Mahajan (2018). Three-dimensional biomechanics of coronary arteries. International Journal of Engineering Science, vol 130, pp. 93-114, https://doi.org/10.1016/j.ijengsci.2018.03.002

2.  S. Mishani, H. Belhoul-Fakir, C. Lagat, S. Jansen, B. Evans, M. Lawrence-Brown, (2021). Stress distribution in the walls of major arteries: implications for atherogenesis. Quantitative Imaging in Medicine and Surgery, vol 11, https://qims.amegroups.com/article/view/68560

3.  G. A. Holzapfel, R. W. Ogden (2018). Biomechanical relevance of the microstructure in artery walls with a focus on passive and active components. American Journal of Physiology-Heart and Circulatory Physiology, vol 315, https://doi.org/10.1152/ajpheart.00117.2018

4.  F. Gervaso, C. Capelli, L. Petrini, S. Lattanzio, L. Di Virgilio, F. Migliavacca (2016). On the effects of different strategies in modelling balloon-expandable stenting by means of finite element method. Journal of Biomechanics, vol 41, pp. 1206-1212, https://doi.org/10.1016/j.jbiomech.2008.01.027

5.  T. Ramanathan, H. Skinner, (2005). Coronary blood flow. Continuing Education in Anaesthesia Critical Care & Pain, vol 5, pp 1743-1816, https://doi.org/10.1093/bjaceaccp/mki012