Example 9a2: Advection-diffusion coupled to gas exchange in asymmetric airways with four lumped acinar units

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 9a2: Calculate advection diffusion gas-exchange solution.
# 
# Required data:
#   Radius field: Either read field file contains airway radii or set up using 
#     "fem update mesh" command.
#   Ventilation field: Read field file containing fractional ventilation values 
#     at every node
#   Perfusion fields: Read field file containing fractional perfusion, capillary
#     blood volume (mm3), RBC transit time (s). Only terminal/acinar nodes req'd.
# 
#   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;
}

#=================================
## Files
$geom="test_geom";
$radius="test_radius"; #field1=radius(mm)
$V_file="test_ventilation"; #field2=a/A ratio, #field3=V(fractional)
$Q_file="test_perfusion"; #field4=Q(fractional), field5=Vc(mm3), field6=tt(sec)
# system "rm test.error test.nodes test.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 = 10;
$T_insp =2.0; #seconds
$T_expn = $T_insp; #seconds
$acinus = 10; #initial acinar volume (mm3/acinus)
$flow = 20; #at inlet (mm^3.s-1)
$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 node;r;$geom;
fem de elem;r;$geom 1d;
fem update scale_factor normalise;
fem define mesh calculate #sets up versions for 1D tree & element fields
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 de field;r;$radius values start $nj_radius end $nj_radius nj_radius $nj_radius;
fem update field $nj_aA substitute constant 1; #aA ratio s/b 0 for all conducting airways
fem de field;r;$V_file values start $nj_V end $nj_V; # read in flows
fem update field $nj_V summation normalise
fem def field;r;$Q_file values start $nj_Q end $nj_tt; 

#=================================
## Update problem
fem de moti;c lung;
fem de;add mesh;c airway lumped_parameter volume $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;
}

Files used by this example are:

Name                     Modified     Size

example_9a2.com          30-Mar-2011  5.0k
GasExchange.ipequa       30-Mar-2011  2.1k
GasExchange.ipsolv       30-Mar-2011  2.0k
GasExchange.ipvalu       30-Mar-2011  436 
Linear_Line.ipbase       30-Mar-2011  893 
N2O2.ipmate              30-Mar-2011  288 
draw.com                 30-Mar-2011  1.6k
test.ippara              30-Mar-2011  5.8k
test_geom.ipelem         30-Mar-2011  2.2k
test_geom.ipnode         30-Mar-2011  2.0k
test_perfusion.ipfiel    30-Mar-2011  1.3k
test_radius.ipfiel       30-Mar-2011  2.2k
test_ventilation.ipfiel  30-Mar-2011  2.0k

Download the entire example:

Name                      Modified     Size

examples_9_9a_9a2.tar.gz  06-Apr-2011  120k

Testing status by version:

Testing status by file:

• status
  

  i686-linux
  1	cm:	error exit status 1.
  1	cm-debug:	error exit status 1.
  x86_64-linux
  1	cm:	error exit status 1.
  1	cm-debug:	error exit status 1.

  
  

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

Input last modified: Wed Apr 6 17:28:06 2011