Example 9a3: Advection-diffusion coupled to gas exchange in whole lung: symmetric airways with lumped acinar unit

This example predicts the oxygen concentration distribution in a simple symmetric (single path) airway model geometry by solving the 1D advection-diffusion equation coupled to gas exchange ODEs in the acinus.

Created by Merryn & Annalisa 06 April 2011.


The comfile run by this example is as follows:

########################################################
# Example 9a3: Calculate advection diffusion gas-exchange in whole lung
#              with 16 generations of symmetric conducting airways.
# 
# General setup.
#   Radius field: Symmetric mesh, so uses Weibel generation data from
#     "fem update mesh" command.
#   Ventilation field: Uniform distribution calculated using 
#     "fem evaluate flow uniform" command.
#   Perfusion fields: Set to be single values at the one acinus.
# 
#   Solve over breath/s (multiple breaths req'd to reach steady-state).
########################################################

# If the example path is not set, default to current directory
if (! defined $example) {
    $example = "./";
}
# Drop off the trailing / in the example path
$chopped = chop $example;
if ($chopped ne "/") {
    $example .= $chopped;
}

#=================================
# system "rm test.error test.nodes test.history"; 
system "mv test.history last.history"; 

#=================================
## Fields (defined in ipequa file)
$nj_radius=1; # set up in fem de mesh; c airway field
$nj_aA=2;     # duct to alveolar cross-sectional area ratio
$nj_V= 3;     # normalised (fractional) ventilation per unit
$nj_Q=4;      # normalised (fractional) perfusion per unit
$nj_Vc=5;     # capillary blood volume
$nj_tt=6;     # RBC transit time
$nj_pob=7;    # capillary blood PO2
$nj_poa=8;    # alveolar air PO2
$nj_source=9; # flux source term
$num_fields=9;

#=================================
## Parameters
$INSP=1;
$EXPN=2;
$FIXED_NODE = 1;
$FLUX_NODE = 1;
$NUM_BREATHS = 20;
$T_insp =2.0; #seconds
$T_expn = $T_insp; #seconds
$acinus = 3e6/(2**15); #acinar volume for a lung with FRC=3 litres
$flow = 2.0e5; #at inlet (mm^3.s-1), tidal volume = 500 ml
$initial = 5.8E-6; #initial concentration ~= 105mmHg (mmol/mm^3)
$conc = 8.26E-6; #inspired concentration ~=150 mmHg (mmol/mm^3)
$dt = 0.01; 

#=================================
## Set up problem
fem de param;r;test;
fem de coor 3,1;
fem de base;r;Linear_Line;

#=================================
## Geometry
fem de mesh;r;SymmetricConducting;
fem de mesh;c airway field num_field $num_fields noversions 5..10;
fem evaluate ordering; 
fem group elem all by terminal as TERMINAL;
fem group node elem TERMINAL xi1=1 as TERMINAL;

#=================================
## Fields
fem update mesh geometry radius_field $nj_radius anat weibel;

#=================================
## Update problem
fem de moti;c lung;
fem de;add mesh;c airway lumped_parameter volume $acinus;
fem evaluate flow uniform; #note that this overwrites $nj_V;
# note that the $nj_aA field MUST be set AFTER evaluate flow, otherwise it will contain junk 
fem update field $nj_aA substitute constant 1; #aA ratio s/b 0 for all conducting airways
fem update field $nj_Q substitute constant 1/(2**15);
#fem update field $nj_Vc substitute constant 6.3;  # capillary blood volume per acinus
#fem update field $nj_tt substitute constant 0.75; #transit time per acinus

##########################################################
# NEED TO BE CAREFUL ABOUT BALANCING THE AIR-BLOOD SURFACE AREA,
# THE CAPILLARY BLOOD VOLUME PER ACINUS, AND THE TRANSIT TIME

# note that the perfusion model gives air-blood surface area
# per acinus as an output, but this isn't read in here. This should be
# included in future.

# FUNCTIONAL STUDIES (WITH DECRECRUITED CAPILLARIES) GIVE ~ 100 ML CAP BLOOD
# MORPHOMETRIC MEASUREMENTS GIVE ~ 200 ML CAP BLOOD

# INCREASING THE CAPILLARY BLOOD VOLUME GIVES SIGNIFICANTLY LOWER PAO2
# INCREASING THE TRANSIT TIME SIGNIFICANTLY INCREASES THE PAO2
# INCREASING THE AIR-BLOOD SURFACE AREA SIGNIFICANTLY INCREASES THE PA02

$capillary_tt = 0.8; #seconds
#$cardiac_output = 83333; #mm^3/sec (== 5 litres/min; set in GasExchange.ipvalu)
$cardiac_output = 1e5; #mm^3/sec (== 6 litres/min; set in GasExchange.ipvalu)
# implies that the blood volume is:
$blood_volume = $cardiac_output * $capillary_tt;

fem update field $nj_Vc substitute constant $blood_volume/(2**15);  # cap blood volume per acinus
fem update field $nj_tt substitute constant $capillary_tt; #transit time per acinus

fem de equa;r;GasExchange class $INSP;
fem de equa;r;GasExchange class $EXPN;
fem de valu;r;GasExchange; # contains parameter values

#=================================
#...DEFINE THE PROBLEM SETUP FOR INSPIRATION
fem de mate;r;N2O2 class $INSP;
fem update motion inlet_flow $flow class $INSP;
fem de init;c lung initial $initial class $INSP;
fem de;add init;c lung fixed node $FIXED_NODE conc $conc class $INSP;
fem de solve;r;GasExchange name test class $INSP;

#...DEFINE THE PROBLEM SETUP FOR EXPIRATION
fem de;add init;c lung flux node $FLUX_NODE class $EXPN; 
fem de;add init;c lung fixed node TERMINAL class $EXPN;
fem de solve;r;GasExchange name test class $EXPN;

#=================================
# #FOR EACH BREATH
$T_start = 0;
$T_end = $T_insp;
$restart="";
for $breath (1..$NUM_BREATHS){

## Inspiration
      if($breath > 1){
            $T_start = $T_end;
            $T_end = $T_end + $T_insp;
            fem update solution class $INSP substitute solution class $EXPN;
            fem update solution node $FIXED_NODE class $INSP substitute constant $conc;
            fem update motion inlet_flow $flow dv_from $EXPN class $INSP;
            fem update mass from $EXPN class $INSP;
      }
      fem solve $restart from $T_start to $T_end delta_t $dt class $INSP;
      if($breath==$NUM_BREATHS){
            fem export node;test_end_insp as end_insp;
            fem export elem;test_end_insp as end_insp;
      }
      fem update solution class $EXPN substitute solution class $INSP;

## Expiration
      $T_start = $T_end;
      $T_end = $T_end + $T_expn;
      $restart="restart";
      fem update motion inlet_flow -$flow dv_from $INSP class $EXPN;
      fem update mass from $INSP class $EXPN;
      fem solve $restart from $T_start to $T_end delta_t $dt class $EXPN;
      if($breath==$NUM_BREATHS){
            fem export node;test_end_exp as end_exp;
            fem export elem;test_end_exp as end_exp;
      }
      fem update solution class $INSP substitute solution class $EXPN;
}

fem quit;

Files used by this example are:

Name                        Modified     Size

example_9a3.com 06-May-2011 6.0k GasExchange.ipequa 06-Apr-2011 2.1k GasExchange.ipsolv 06-Apr-2011 2.0k GasExchange.ipvalu 06-May-2011 423 Linear_Line.ipbase 06-Apr-2011 893 N2O2.ipmate 06-Apr-2011 288 SymmetricConducting.ipmesh 06-May-2011 2.4k draw.com 06-Apr-2011 1.6k test.error 06-May-2011 117k test.history 06-May-2011 320k test.ippara 06-Apr-2011 5.9k test.nodes 06-May-2011 398k test_end_exp.exelem 06-May-2011 65k test_end_exp.exnode 06-May-2011 195k test_end_insp.exelem 06-May-2011 65k test_end_insp.exnode 06-May-2011 195k

Download the entire example:

Name                      Modified     Size

examples_9_9a_9a3.tar.gz 07-May-2011 145k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmFailureSun Mar 6 00:05:15 20160
last breakSun Apr 17 00:48:00 20111
cm-debugFailureSun Mar 6 00:05:16 20161
last breakThu Apr 7 00:07:00 20111
x86_64-linux
cmFailureSun Mar 6 00:01:26 20161
last breakSun Apr 10 00:30:00 20112
cm-debugFailureSun Mar 6 00:01:26 20161
last breakThu Apr 7 00:05:00 20110

Testing status by file:


Html last generated: Sun Mar 6 05:50:37 2016

Input last modified: Fri May 6 10:17:55 2011


CMISS Help / Examples / 9 / 9a / 9a3