Frog Ventricular Pressure-Volume Curve Model
Contents
Description
- In this exercise you will set up and run a nonlinear finite elastic analysis of the passive inflation of a truncated ellipsoidal model of a frog left ventricle.
- You will need measurements of the long axis (base-apex) outer dimension, short axis (equatorial) outer diameter, and equatorial wall thickness all in mm to customize your model geometry.
- The mesh consists of three elements that are wrapped onto themselves in a prolate spheroidal coordinate system. The basis functions are linear in the circumferential direction and cubic in the longitudinal and radial directions (linear-cubic-cubic). The material chosen is a transversely isotropic exponential strain energy function. The material is stiffer in the fiber direction than perpendicular to them. The fibers whose orientation is defined by fiber angles in the node form with respect to the circumferential direction vary linearly in radial direction across the wall from -37 degrees on the epicardium to +84 degrees on the endocardium. The boundary conditions assure that all rigid body motions are suppressed. Pressure is prescribed on the endocardium to simulate the passive filling of a ventricle (diastole).
- Model output will include a table of ventricular pressures and volumes, which you will try to match to experimental measurements by adjusting one of the material coefficients (stress scaling coefficient) of the constitutive model.
The Continuity model file provided already has an approximate geometry, fiber angles and initial material properties prescribed. It has also been set up with an initial pressure increment per time step of 0.5 mm Hg. If you have to change the geometry or material properties substantially, it is possible that the model may fail to converge especially at low pressures. In that case you may need to lower the pressure increment in the Systemic Circulation tab of the Circulation Model.
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
Load and Render the Initial Model
File→Load Model... or click
Select the cont6 file for this tutorial (download prolatefrogPV.cont6 )
Click Open.
Mesh→Calculate Mesh... or click
Click OK to close the Calculate Mesh Form and execute mesh calculations on the server
Mesh→Render→Elements... or click
Click the lines radio button
Click Render to display mesh lines
- The mesh should now look similar to the screenshot below. If you like you can also render fibers and surfaces of this model.
Customize Ventricular Geometry Using Measurements
You will need your measurements in mm of the outer (epicardial) apex-base length (L), outer (epicardial) short-axis (equatorial) diameter (D) and equatorial wall thickness (H). From these you will calculate: Prolate Spheroidal Focus Position (d) in mm; inner (endocardial) Lambda coordinate; and outer (epicardial) Lambda coordinate as shown below (note that asinh below is the arcSinh function):
First customize the prolate spheroidal model by specifying the focus position, d that you calculated
Note that prolate spheroidal is already selected in the Global Coordinates: pop-up menu
Replace the focal length currently in the focus position input field with the value of d that you calculated for your frog heart.
Click OK to submit Coordinate Form
- Next specify the inner and outer Lambda coordinates of your prolate spheroidal model using the values of lambda that you calculated:
Note in Node List that nodes 1, 2, 5 and 7 have the name ENDO in their label shown at the top of the form. These are the endocardial (inner surface) nodes. Select each one of these nodes in turn, and under Coordinate 1 in the field next to Value, replace the default entry with the value you calculated for Lambdaendo. Press return after entering the new value.
Note in Node List that nodes 3, 4, 6 and 8 have the name EPI in their label shown at the top of the form. These are the epicardial (outer surface) nodes. Select each one of these nodes in turn, and under Coordinate 1 in the field next to Value, replace the default entry with the value you calculated for Lambdaepi. Press return after entering the new value.
Click OK to submit Node Form
- Next you need to specify a point on the plane of the base to ensure that cavity volumes are calculated correctly:
Under Define the Origin replace the value next to X with -L/3. If you measured L as 6 mm, then the current value of -2.0 will be correct.
Click OK to submit List Volume Form and you will see the volume of the ventricular cavity and wall listed
- It is possible that your unloaded volume may differ between the model and experiment if the unloaded volume is not consistent with your geometric measurements, but this is more likely to reflect an inaccuracy in your volume measurement that your geometry.
- Be sure to update the initial conditions and save:
Biomechanics→Update→Initial conditions with undeformed nodes
This command copies the nodal geometric coordinates and derivatives from the Node Form to the Initial Conditions tab of the Boundary Conditions Form
Mesh→Calculate Mesh... or click
Click OK to close the Calculate Mesh Form and execute mesh calculations on the server
File→Save→Save model as... and save your model as a .cont6 file with a new name, in case you need to restart later.
- Render the new geometry and clean up:
Mesh→Render→Elements... or click
Click the lines radio button
Click Render to display your new mesh
View→Show→OpenMesh... or click on
Click on 1. element lines1 in the list on the left, and click the Delete button beneath the list to remove the rendering of the old mesh.
- Close the window.
Run the Biomechanics Model
Biomechanics→Solve Nonlinear...
Set Number of Steps to twice the maximum ventricular pressure you reached in your experiment in mm Hg. i.e., If you loaded your ventricle to 25 mm Hg, set Number of Steps to 50, then click the Solve button, and wait for the solver to finish. It may take several minutes. While the solution is being computed, the window will remain open. There will also be output listed to the console window and the Python shell which can be viewed by clicking on
This will initially attempt to solve the model in pressure steps of 0.5 mm Hg. If the solve succeeds, you will see a message in the Solve Nonlinear Form message area. The message will include text saying: The output files from this solve can be found in /Users/Yourname/.continuity/working/ (for Mac OS X) or The output files from this solve can be found in C:\Documents and Settings\Your Directory\.continuity/working/ for Windows.
Locate the file named hemodata.xls in this directory, move it to a more convenient directory and open it with Excel. Column D will contain ventricular pressures in mm Hg and Column E will contain ventricular volumes in mm3. Use Excel to plot the pressure-volume curve and compare it with your measurements.
- If the solution fails to converge or produces an error it may be that you need to reduce the pressure increment:
Biomechanics→Edit→Circulation Model...
Click Ok to edit the current model
Select the Systemic Circulation tab and enter a smaller value for the Pressure increment per time step. This will mean that you will need more time steps to reach the maximum ventricular pressure of your experiment.
File→Save→Save model to save your model, restart Continuity, reload the model, calculate mesh, and try again.
Adjust the Material Properties
Start by resetting Continuity (File→Reset; click on all reset and all lose data options), then re-load your previous model file.
Biomechanics→Edit→Constitutive Model...
Click on the Set parameters tab, select stress_scaling_coefficient from the list on the left and change the entry after Value in the right panel. Increasing this coefficient will increase the stiffness of the myocardium and proportionately increase the slope of the P-V curve. So for example, if your model predicted pressures that were about 50% of expected values for given volumes, then double the value of the coefficient. If you have to reduce the coefficient significantly, you will probably need to decrease the pressure step as described above for the solution to converge.
Click on the Edit Equations tab and then click OK
Mesh→Calculate Mesh... or click
Click OK to close the Calculate Mesh Form and execute mesh calculations on the server.
Biomechanics→Solve Nonlinear...
Set the appropriate number of time steps, click the Solve button, and wait for the solver to finish again.
- Compare the model P-V curves with your experimental data and if necessary adjust the material parameter and solve again.
Render and List Stress and Strain Solutions
Biomechanics→Render→Surface...
Change the pop-up menu choice after At Xi to 1 and change Location to 0.0
Select the deformed radio button
Under Variables select E-out- Lagrangian Green's Strain Tensor wrt material coordinates
Click OK to render the fiber strain
Biomechanics→Render→Surface...
Change the pop-up menu choice after At Xi to 1 and change Location to 0.5
Select the deformed radio button
Under Variables select E-out- Lagrangian Green's Strain Tensor wrt material coordinates
On the components matrix, click ZZ (ZZ is [2, 2] in the newest version of Continuity)
Click OK to render the radial strain
- The result should look something like the screenshot below:
For BENG 172, you need to include a similar screenshot that shows the fiber and radial stresses vs. radius instead of the strains vs. radius. Remove the renderings of the strain distributions, and repeat the above steps. This time select stress-out - Second Piola-Kirchhoff Stress Tensor wrt material coordinates. Please include a scale bar along with your screenshot.
View→Show→OpenMesh... or click on
In the Colors tab, you will see the range that the color scheme varies between (scale bar).
Biomechanics→List Stress and Strain...
- Note: This step is optional for BENG 172.
In the Variables tab, use checkboxes to deselect any variables you do not wish to list
Note that the Output Variables are calculated using the equations entered in the Output Variable folder of the submitted Constitutive model
In the Locations tab, you can selected the element list and points which which solutions will be listed
Click OK to display a listing of the selected Output Variables in the Table Manager
Make a movie
- Note: This step is optional for BENG 172.
- Once you have achieved a satisfactory match between the model and your experimental measurements, you can make an animation of your model filling, and render or list stress and strain distributions. Follow these steps to make the movie:
Save you model and restart Continuity
Mesh→Calculate Mesh... or click
Click OK to close the Calculate Mesh Form and execute mesh calculations on the server.
Mesh→Render→Elements... or click
Click the lines radio button
Click Render to display undeformed mesh
Biomechanics→Solve Nonlinear...
Leave Number of Steps as 1 and click Solve.
After the solution has finished, leave the Solve Nonlinear Form open.
Mesh→Render→Elements... or click
Click the deformed radio button
Click Render to display deformed mesh lines
Back in the Solve Nonlinear window click Solve again, and after it has finished, again render deformed lines.
Now increase the Number of Steps in the Solve Nonlinear Form to 2 and repeat a couple more solves and renders
Now increase the Number of Steps in the Solve Nonlinear Form to 4 and repeat solving and rendering until you reach maximum pressure.
- View→Animate... to open the Animation Player
Click the SetAnim button to open the SetAnimation window
Click the AddFrame button once for each time you rendered a solution plus one more for the undeformed state (Frame 0)
Double click on Root on the list on the left to expose the list of element line objects that you rendered
One at a time, select each frame from the list on the right, select the corresponding element line object starting with the undeformed mesh ~element lines1 and click the AddGeom button. Repeat adding each successive element lines object to a new frame until you have assigned each lines object to a separate frame.
Click OK to close the SetAnimation window.
Click the play button
to play back the animation. Try rotating the model while the playback loops.