Inflation of a Transversely Isotropic Thick-Walled Cylindrical Tube

Contents

  1. Inflation of a Transversely Isotropic Thick-Walled Cylindrical Tube
    1. Description
    2. Start Continuity
    3. Create Mesh
    4. Refine the Mesh
    5. Add Biomechanics Data
    6. Solve Biomechanics Problem
    7. Calculate and Render Strains
    8. Pre-built model

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

  • Mesh→Edit→Coordinates…

    • Select cylindrical polar in the Global Coordinates: pop-up menu

    • Click 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 compilation

    • Click OK to submit the Material Coordinate model

  • Mesh→Edit→Basis…

    • 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

  • Mesh→Edit→Nodes…

    • Click Insert Node in the left panel 3 times to create a total of 4 nodes

    • Note 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

  • Mesh→Edit→Elements…

    • 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, 4in 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 button

    • Click Render to display the mesh as a wireframe

  • Mesh→Render→Elements… or click 

    • Click the surfaces radio button

    • Click 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.

  • 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.

  • File→Library→Search…

    • 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 OKbutton 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…

    • Mesh→Edit→Nodes…

    • After having opened the two forms, confirm that the Boundary Conditions forms has the same basis functions as the Nodes form

    • Close the Nodes form

  • Biomechanics→Edit→Boundary Conditions…

    • Click on the Deformed Coordinate 2 tab

      • Click the Insert Nodes button

      • Enter .*proximal under Nodes(s): and 0.0 under Value:

      • Click the Insert Nodes button

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

    • Click on the Deformed Coordinate 3 tab

      • Click the Insert Nodes button

      • Enter .*proximal under Nodes(s): and 0.0 under Value:

      • Click the Insert Nodes button

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

    • Click on the External Pressure tab

      • Change 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 job

    • Note 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

  • Mesh→Render→Elements…

    • Click the lines radio button

    • This time select the deformed radio button

    • Click 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,Bentry 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 element

    • Click the surfaces radio button

    • Click the deformed radio button

    • Click 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 button

    • under 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

    • View→Show Open Mesh …

    • 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 T

    • In 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 solution

    • Click 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

  • File→Library→Search…

  • 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.