## Inflation of a Transversely Isotropic Thick-Walled Cylindrical Tube

Contents

### Description

- This example guides you through setting up and running the geometrically and materially nonlinear problem of the inflation of a thick-walled cylindrical tube with orthotropic exponential material properties.
- The constitutive model is obtained from published experimental data (doi:10.1152/ajpheart.01323.2005).
- The basic mesh consists of one cylindrical element wrapped onto itself defined in the cylindrical polar coordinate system. The basis functions for deformed and undeformed coordinates are trilinear Lagrange with 27 Gauss quadrature integration points. This mesh is refined to contain ten trilinear Lagrange elements.
- Inflation is achieved by applying a pressure on the inner surface of the mesh. The other boundary conditions assure that rigid body motions are suppressed.
- Even using ten trilinear Lagrange elements, we can not assume the solution is well converged. The solution in the center of the elements will be most accurate. As we refine the mesh and get more elements and/or use higher-order interpolation, the solution will probably change as it approaches the exact solution.

### Start Continuity

- Launch the Continuity 6.4 Client
On the About Continuity 6.4 startup screen

Leave the

**mesh**checkbox checked under`Use Modules:`In addition, check the

**biomechanics**checkbox

Click

**OK**to bring up the main window

### Create Mesh

Select

**cylindrical polar**in the`Global Coordinates:`pop-up menuClick

**OK**to submit`Coordinate Form`

Mesh→Edit→Material Coordinates...

Use the equation below to write

**dY_dMatl**, the material coordinate transformation from cylindrical polar to rectangular cartesian.- dY_dMatl = Matrix([[cos(X[1]), -sin(X[1]),0],[sin(X[1]),cos(X[1]),0],[0,0,1]])
The SymPy tutorial demonstrates how to create a matrix in SymPy

Select

**Compile...**and click**Compile Code as: C Double Precision**Click

**Ok**to close the notice for a successful compilationClick

**OK**to submit the`Material Coordinate`model

Choose

**Lagrange Basis Function→3D→Linear-Linear-Linear**Click

**Add Linear-Linear-Linear**- Verify that the list of basis functions now contains:
- Linear-Linear-Linear Lagrange 3*3*3
- Linear-Linear Lagrange 3*3

Click

**OK**to submit`Basis Form`

Click

**Insert Node**in the left panel 3 times to create a total of 4 nodesNote that the first basis function you created (

**Linear-Linear-Linear Lagrange 3*3*3**) with the`Edit Basis Function`command is already selected in the menus below`Coordinate 1`,`Coordinate 2`, and`Coordinate 3`In the

**Value**fields under`Coordinate 1`,`Coordinate 2`, and`Coordinate 3`enter the following (R,Theta,Z) nodal coordinates:- Node 1: (1., 0., 0.)
- Node 2: (1., 0., 4.)
- Node 3: (1.5, 0., 0.)
- Node 4: (1.5., 0., 4.)

Replace the

`default`node label with the appropriate labels below.**Note:**You must hit enter for the label to be updated.- Node 1: inner_proximal
- Node 2: inner_distal
- Node 3: outer_proximal
- Node 4: outer_distal

Enable a

`Field Variable`by selecting the**Field Vector 1**tab and choosing**Linear-Linear-Linear Lagrange 3*3*3**from the`Select Basis Number`menu below`Field variable 1`.Click

**OK**to submit`Nodes Form`

Element 1 has eight element nodes, but only four global nodes because it wraps around on itself in the theta direction, which corresponds to

`xi1`. Enter**1, 1, 2, 2, 3, 3, 4, 4**in the`Global Node Numbers`boxes. Use the`tab`key to change the input focus to the next box. Note that the order that global node numbers are entered determines the local`Xi`coordinate directions in the element, as illustrated by the graphic in the input form.Click

**OK**to submit`Element Form`

File→Send or click

Mesh→Calculate Mesh... or click

Mesh→Render→Nodes... or click

Mesh→Render→Elements... or click

Click the

**lines**radio buttonClick

**Render**to display the mesh as a wireframe

Mesh→Render→Elements... or click

Click the

**surfaces**radio buttonClick

**Render**to display inner mesh surfaces (at`Xi(3)`=0.)- The mesh should now look similar to the screenshot below

### Refine the Mesh

To get a sufficiently converged result using linear elements, it is necessary to use multiple elements. Therefore, we will refine our single element into many elements.

Decrease the

`New Element per old element in``xi1`and`xi2`to 1. Increase the`New Element per old element in``xi3`to 10Click

**OK**to refine the mesh

File→Send or click

Mesh→Calculate Mesh... or click

Mesh→Render→Nodes... or click

Mesh→Render→Elements... or click

Click the

**lines**radio buttonClick

**Render**to display the mesh as a wireframe

To delete or hide previously rendered objects View→Show→OpenMesh... or click on

Note that the element labels are preserved (i.e. all elements at R=0 have

`proximal`in their label)

### Add Biomechanics Data

A biomechanics problem requires material properties and boundary conditions before it can be solved. A nonlinear hyperelastic constitutive equation from literature will be used for the material properties. The boundary conditions will define how the biological testing is preformed.

##### Incorporate a Published Constitutive Model

An orthotropic porcine (pig) coronary artery model will be used (Wang et. al. Three-dimensional mechanical properties of porcine coronary arteries: a validated two-layer model. *Am J Physiol Heart Circ Physiol* 291:H1200-H1209, 2006 doi:10.1152/ajpheart.01323.2005). An exponential constitutive model template from the continuity database will serve as the basis for the model. While the experimental paper presumes incompressibility, the template allows for slight incompressibility governed by a bulk modulus. This improves numerical stability.

*Load a constitutive model from the Continuity database*Select model 1097:

**BM_InflateTube_ConstitutiveModel_Basic**Right click on the model title and select

**Load**In the

`Resolve Conflict`window that appears, select the radial button for**Retain current problem, but overwrite the following object:**Click

**OK**to load this constitutive model from the database

Biomechanics→Edit→Constitutive Model...

*Edit the constitutive model*All of the strain calculations (detailed in the StrainIn3DBlock tutorial) are prepopulated in this template. The stress calculations, including the compressibility are, are already included as well. The only variable that needs editing is

**Q**.This is a exponential (Fung type) constitutive equation and

**Q**(under`Calculated Variables→energy`) is the exponential term. Equation A6 in the appendix of the reference will be used for**Q**.**Q = <b1>*E[1,1]*E[1,1] + <b2>*E[2,2]*E[2,2] + <b3>*E[0,0]*E[0,0] + <b5>*2.0*E[2,2]*E[0,0] + <b6>*2.0*E[1,1]*E[0,0] + <b4>*2.0*E[1,1]*E[2,2]**

The stress calculation combines two strain energy terms. While the first term,

**W1**is well known, the second ,**W2**, relaxes the incompressible assumption and allows the**bulk modulus**to dictate the change in volume. This makes the problem easier to solve numerically.Note that the coefficients defined in <angled brackets> appear in the

`Set parameters`tab and can be changed without re-compiling. Additionally, Python indexes from zero, so the first component of the strain tensor is E[0,0] and not E[1,1]. Using the above material coordinate transformation, the ordering of the strain tensor components is r, theta, z corresponding to 0, 1, 2 in python indexing.Select

**Compile...**and click**Compile Code as: C Double Precision**- The following output variables need to be added.
- Right-click on the variable stress_out in the left panel of the Edit Equations tab and select Insert variable here...
Set the variable type to

*Evaluated*- T = (F*stress*F.T)/J - Cauchy Stress Tensor
- N = stress*F.T - Nominal Stress Tensor

Click to the

**Set parameters**tab to change the`Value`for each coefficient**b1**- 1.25**b2**- 2.46**b3**- 0.55**b5**- 0.06**b6**- 0.08**b4**- 0.36**stress_scaling_coefficient**- 4.46**bulk_modulus**- 350

Click

**OK**to clear the notice for a successful compilation and select the**OK**button to submit the constitutive model

#### Define the Boundary Conditions

The boundary conditions for this model will apply a pressure on the inner surface and fix the open ends of the vessel causing the inflation of the tube while preventing rigid body motion.

Biomechanics→Update→Initial conditions with undeformed nodes

Biomechanics→Edit→Boundary Conditions...

After having opened the two forms, confirm that the

`Boundary Conditions`forms has the same basis functions as the`Nodes`formClose the

`Nodes`form

Biomechanics→Edit→Boundary Conditions...

Click on the

**Deformed Coordinate 2**tabClick the

**Insert Nodes**buttonEnter

**.*proximal**under`Nodes(s):`and**0.0**under`Value:`Click the

**Insert Nodes**buttonEnter

**.*distal**under`Nodes(s):`and**0.0**under`Value:`

Click on the

**Deformed Coordinate 3**tabClick the

**Insert Nodes**buttonEnter

**.*proximal**under`Nodes(s):`and**0.0**under`Value:`Click the

**Insert Nodes**buttonEnter

**.*distal**under`Nodes(s):`and**0.0**under`Value:`

Click on the

**External Pressure**tabChange

`Select Pressure Type`to**Incremental**Specify the pressure increment to be

**0.4**on the**Inner Surface**- The units of the pressure increment depend on your mesh scale and constitutive equation parameters.

- It is a good idea to now go back through the Boundary Conditions Form to double check the parameters you just set up
Click

**OK**

### Solve Biomechanics Problem

File→Send or click

Mesh→Calculate Mesh... or click

Click

**OK**to close the`Calculate Mesh Form`and execute mesh calculations on the server

Biomechanics→Solve Nonlinear...

For

`Time Step`, enter**0.1**For

`Number of Steps`, enter**10**Click the

**Solve**button, and wait for the solver to finish its jobNote that the time counter (

`Initial time`) updates to 1.0 after the solve. For ever 1.0 increase in the time counter, the boundary conditions are implemented once. For example, the pressure is incrementally increased to 0.4 during this simulation. If a displacement boundary condition of 0.1 had been given and the simulation run to a total time of 3.0 (`Time Step`*`Number of Steps`), a total displacement of 0.3 would be applied.

### Calculate and Render Strains

Click the

**lines**radio buttonThis time select the

**deformed**radio buttonClick

**Render**to display deformed mesh lines

View→Show→OpenMesh... or click on

Click on

**3. element lines2**in the list on the left, and enter**1,0,0**in the`R,G,B`entry field to change the mesh lines from blue to red.- Press [return] and close the window
- The mesh should look like the screenshot below

Mesh→Render→Elements... or click

In the

`Element List`, enter**1**to render the inner surface of the elementClick the

**surfaces**radio buttonClick the

**deformed**radio buttonClick

**Render**to display inner mesh surfaces (at`Xi(3)`=0.)

Biomechanics→Render→Surface...

Change the pop-up menu choice after

`At Xi`to**2**and change`Location`to**0.5**Select the

**deformed**radio buttonunder

`Variables`select**T Cauchy Stress Tensor**Click OK

**to render the fiber (circumferential) strain**

- The mesh should now look similar to the screenshot below

- If you do not get the exact color- try doing this
Select Textured Field4

**from the**Rendering order**list**Click on Colors

**tab select**RGB**or**BGR**for**Color Scheme:- Close the window

Biomechanics→List→Stress and Strain...

Unselect All Variables

**and reselect**X**and**TIn the Locations

**tab, select**Number of Points**and change**`xi1`,`xi2`and`xi3`to 1. This will export the solutions at the center of each element (`xi1`,`xi2`and `xi3' = 0.5) and provide the radial stress distribution. Since this model uses only a small number of linear finite elements, the center of each element has the most accurate solutionClick OK

**to display a listing of the selected**`Output Variables`in the Table Manager- This table can be saved (File→Save...) and graphed using external software
- Please note that the solution is most accurate at the middle Gauss point within each element particularly with using linear-linear-linear basis functions

- Biomechanics→List→Nodal Solutions...
- The nodal solutions table contains the final node locations and the residual values
- Since Continuity minimizes these residuals when solving the equilibrium equations, they will be very small. However, if a displacement boundary condition was enforced, the residual will be non-zero. The residual will be the force experienced by the body to achieve that displacement
- In this model, the residuals in the Z direction are non-zero. This is a result of the zero displacement boundary conditions in the Z direction
- The opposite of the residuals could be applied as a force boundary condition to achieve the same displacement

### Pre-built model

- A pre-built model for this tutorial is available on Continuity's public database.
Title:

**bm_inflate_tube**ID:

**1387**

A script is provided in the script manager

*to run the tutorial**Note that the material coordinate model and constitutive model may still need to be compiled.*