A single prolate. Geometry is interpolated isochorically
(Focus^3.Cosh[lambda].Sinh^2[lambda]). Hydrostatic pressure is interpolated
quadratically in Xi3 (element based) and the pole-zero law defines the
incompressible material behaviour. The prolate is inflated then contracted isovolumically then ejection is simulated.
#Example_5532 Inflation and contraction of an
# orthotropic prolate spheroid with fibres and sheets
fem
fem define coordinate;r;profib;example #Defines 3 GLOBAL PROLATE SPHEROIDAL
# ! coordinates.
# ! Transmural interpolation is in
# ! Focus^3.Cosh(Lamda).Sinh^2(Lamda).
fem define node;r;;example #Focus = 37.5. Reads 4 NODES, 3 COORDS,
# ! 1 VERSION per node per nj, 0
# ! DERIVATIVES for all coords. Nodal
# ! coordinates, Xj, are (lamda,mu,theta):
# ! 1=(0.43,0,0);2=(0.43,120,0);
# ! 3=(0.72,0,0); 4=(0.72,120,0).
fem define base;r;;example #Defines 3 basis functions:
# ! 1) Lagrange/Hermite with 3 Xi coords,
# ! Linear Lagrange interpolant and 1
# ! Gauss point in Xi1 dirctn (theta-
# ! circumferential) but 3 in Xi2 (mu-
# ! longitudinal) and Xi3 (lamda-
# ! transmural);
# ! 2) An auxiliary basis with 5 auxiliary
# ! element parameters, pressure basis
# ! same number of Gauss points as for
# ! basis function 1, and polynomial
# ! degrees of: 0 for Xi1 and Xi2 for
# ! all 5 parameters; 0 in Xi3 for
# ! parameter 1; 1 for parameter 2; 2
# ! for parameter 3; -1 for parameter
# ! 4 (Xi3=0 face pressure bc); & -2
# ! for parameter 5 (Xi3=1 face
# ! pressure bc).
# ! 3) 2 Xi coords with linear Lagrange
# ! interpolant, and 3 Gauss points
# ! in Xi1 and Xi2 directions.
fem define element;r;;example #Defines 1 prolate spheroidal element.
fem define fibre;r;;example #Defines fibre angle orientation wrt Xi1
# ! direction (theta) in degrees. Uses
# ! basis 1 for fibre angle interpolation
# ! in all elements. No version prompts
# ! or derivs. Nodes 1,2 (inner nodes)
# ! are oriented at 65 degrees while the
# ! outer nodes 3,4 are -55 degrees from
# ! the axis. Defines sheet angles of -45 degrees
# ! at inner nodes 1,2 and +45
# ! degrees at outer nodes 3,4.
fem define element;r;;example fibre #Defines fibre elements.
fem define window #Defines window dimensions in X,Y,Z
# ! directions as (min,max):(-5,5),(-5,5),
# ! (-5,5).
fem draw lines #Makes line segments visible on window.
fem define equation;r;;example #Static 3D finite elasticity problem with
# ! geometric coords as the dependent
# ! variables. Solution is by Galerkin
# ! finite elements in a nonlinear
# ! analysis of a nonlinear equation.
# ! Incompressible material with basis 1
# ! in all elements for dep vars 1,2,3
# ! and basis 2 in all elements for dep
# ! var 4 (hydrostatic pressure).
fem define material;r;;example #Stresses in constitutive law are
# ! referred to body coords. Defines a
# ! hyperelastic constitutive law for
# ! which the strain energy function, W,
# ! is a function of fibre/transverse
# ! strains that is pole-zero in strains
# ! with 16 constant constitutive law
# ! parameters. The Lagrangian strain
# ! components are referred to the
# ! undeformed fibre coordinate system
# ! where the 1-direction is aligned with
# ! fibre axis, 2-direction is the sheet
# ! axis (in the plane of the sheet,
# ! perpendicular to the fibres), and the
# ! 3-direction is the sheet-normal coordinate.
fem define initial;r;;example #Boundary pressure increments are entered
# ! and hydrostatic pressure is matched
# ! across boundaries with adjacent Xi3
# ! faces. Initial displacements are all
# ! zero. For dependent variable 2
# ! (mu-dirn) fix all nodes (1..4);
# ! Dependent variable 3 (theta-dirn),
# ! fix node 2. No force bcs.
# ! Dependent variable 4 (hyd press) fix
# ! element 1, parameter 4 (Xi3=0 face
# ! inside prolate) with increment
# ! magnitude 2, and parameter 5 (Xi3=1
# ! face outside prolate) with increment
# ! magnitude 0. No force bcs.
fem export nodes;heart1 field as heart #export the undeformed nodes
fem export elements;heart1 field as heart #export the undeformed elements
fem define solve;r;profib;example # Read in solution information.
fem solve increment 0.2 iter 3 #Solve the problem.
# !Incrementing of the ventricular
# ! pressure is neccessary to avoid
# ! unstable solutions. This depends on
# ! the constitutive law and pressure
# ! the constitutive law and pressure
# ! magnitude, as the wall is gradually
# ! stiffened.
fem export nodes;heart2 field as heart #export the deformed nodes
fem export elements;heart2 field as heart #export the deformed elements
#
fem solve increment 0.3 iter 3
fem export nodes;heart3 field as heart
fem export elements;heart3 field as heart
#
fem solve increment 0.5
fem export nodes;heart4 field as heart
fem export elements;heart4 field as heart
#
#
#
fem draw lines def dotted #Draw deformed mesh on window.
fem list strain at 1 ref #List strain information wrt reference
# ! coordinates.
fem list stress at 1 ref #List stress information wrt reference
# ! coordinates.!Example_5521 Contraction of an inflated orthotropic prolate element
fem def initial;w;profib_def
fem eval reac
#
#
# Contraction
#
#
# ############################################
# Inflated solution is converged OK!
# ############################################
# Command updated by fixcom.sh on Wed Aug 23 17:49:03 NZT 2000
# Old command: assign variable WALL 1
$WALL = 1
# Command updated by fixcom.sh on Wed Aug 23 17:49:03 NZT 2000
# Old command: assign variable CAVITY 2
$CAVITY = 2
#
# Set up cavity region geometry
fem define region;r;coupled;example
fem define coord;r;cavity;example region $CAVITY
fem define;add base;r;cavity;example
fem define node;r;cavity;example region $CAVITY
fem define element;r;cavity;example region $CAVITY
#
# Set up equations, material properties and initial conditions
# for cavity region
fem define equation;r;coupled;example region $WALL,$CAVITY lock
fem define material;r;cavity_isovol;example region $CAVITY
fem define initial;r;cavity;example region $CAVITY
#
# Define cavity-to-wall coupling and solution procedure information
fem define coupling;r;coupled;example
fem define solve;r;coupled;example coupled region $WALL,$CAVITY
#
# Transfer wall deformation from the inflation problem to
# geometric dependent variables for the cavity region
fem update solution coupled source_region $WALL
#
# Set up the reference state for the cavity region to be the deformed
# configuration from the inflation problem. To set up the
# 'inflated cavity', the x-coord ('in 1' represents the 1st RC coord)
# for the central cavity node (node 6) is set to the average of the
# x-coords for the adjacent interface nodes (node 2 - this could
# be a list of nodes)
fem update solution cavity_reference average 6 in 1 node 2 region $CAVITY
#
# List undeformed (pre-inflation) and inflated cavity volumes
fem list element total region $CAVITY
fem list element deformed total region $CAVITY
# ############################################
# cavity volume inflated from 40ml to 65ml.
# ############################################
#
fem draw lines region $WALL,$CAVITY
fem draw lines dot def region $WALL,$CAVITY
#
# Check that the current coupled solution is converged (as expected!)
fem solve coupled increment 0.0
# ############################################
# inflated state is a converged solution for the coupled problem!
# ############################################
#
# ISOVOLUMIC CONTRACTION
#
# Increment intracell calcium conc and solve for isovolumic contraction
fem define material;r;profib_active;example
fem define active;r;profib_contract;example
fem solve coupled increment 0.0
# ############################################
# The solution procedure hasn't theoretically converged, but
# if you follow the progress of the solver and the final
# constr vs unconstr resids, it's very close to converged.
# Note that the ratio monotonically decreased,
# which is also a sign of a stable solution process.
# ############################################
#fem define initial;r;profib_endisovolcontr;example region WALL,CAVITY
#
# Save deformed configurations for end ejection
fem define initial;w;profib_endisovolcontr region $WALL,$CAVITY
#
fem draw lines deformed dotted rgb=blue region $WALL,$CAVITY
#
# Check that cavity volume has not changed
fem list element deformed total region $CAVITY
# ############################################
# This contraction phase was indeed isovolumic!
# (cavity volume still 65ml).
# ############################################
#
# EJECTION
#
# Decrease cavity 'stiffness parameter' and solve for ejection
fem define material;r;cavity_eject;example region $CAVITY
fem solve coupled increment 0.0
# ############################################
# Ejection problem is SLOW to converge, but note that the ratio
# monotonically decreases, indicating a stable solution.
# May need to repeat the last solve command to get better convergence,
# but will probably not get down to 10^(-10) tolerance.
# ############################################
#
# Save deformed configurations for end ejection
fem define initial;w;profib_endeject region $WALL,$CAVITY
#
fem draw;add lines deformed dotted region $WALL,$CAVITY
#
# Determine cavity volume at end-ejection.
# NOTE: to calc cavity volume, the deformed position of the central
# cavity node must be recalc'ed from the averaged adjacent wall nodes.
# After doing this the current solution is meaningless, ie. best to save
# deformed configurations for both regions before this (as above).
fem update solution cavity_reference average 6 in 1 node 2 region $CAVITY
fem list element deformed total region $CAVITY
# ############################################
# Not a huge ejection fraction (!!!), but neverless a decrease
# in volume from 65ml to 63ml.
# ############################################