Example 9a1: Advection diffusion equation in asymmetric airways with lumped acinar units

This example predicts the gas concentration distribution in a simple asymmetric airway model geometry by solving the 1D advection-diffusion equation. Acini are modelled as lumped volumes at the end of each terminal bronchiole.

Created by Annalisa Swan 18 Feb 2011.

Figure displays the gas concentration distribution with spheres representing acinar units.

To create an image in order to visualise the results of the simulation, use CMGUI with the file draw.com.


The comfile run by this example is as follows:

# Example 9a1
#  Solve advection diffusion over breath/s. Read in flow field file.
#  Asymmetric airways coupled to lumped acinar units at each terminal.

# 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;
}


#=================================
## Fields
$nj_radius=1;
$nj_aA=2;
$nj_flow=3;
$nj_pressure=4;

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

#=================================
## Geometry
fem de node;r;test_geom;
fem de elem;r;test_geom 1d;
fem update scale_factor normalise;
fem define mesh calculate #sets up versions for 1D tree & element field stuff
fem de mesh;c airway field num_field 6;# radius 1;
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 6 substitute constant 0;
fem update field $nj_aA substitute constant 1; #aA ratio s/b 0 for all conducting airways
fem de field;r;ventilation values start 1 end 4; #6; ## read in flows
fem update field $nj_flow summation normalise;


#=================================
## Parameters
$INSP=1;
$EXPN=2;
$FIXED_NODE = 1;
$FLUX_NODE = 1;
$T_insp = 1.0; #seconds
$T_expn = $T_insp;#1.0; #seconds
$acinus_volume = 100; #mm3/acinus;
$flow = 0.3; # mm^3.s-1
$flow_in  = $flow; # mm^3.s-1
$conc = 1.0; # inspired fractional concentration
$initial = 0.5;# initial concentration
$dt = 0.01; 


#=================================
## Update problem
fem de moti;c lung flow_field $nj_flow pressure_field $nj_pressure; #pressure overwritten later
fem de;add mesh;c airway lumped_parameter volume $acinus_volume;
fem de equa;r;AdvectionDiffusion class $INSP;
fem de equa;r;AdvectionDiffusion class $EXPN;


#=================================
#...DEFINE THE PROBLEM SETUP FOR INSPIRATION
fem de mate;r;"N2O2" class $INSP;
fem update motion inlet_flow $flow_in 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;AdvectionDiffusion name AdvectionDiffusion 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;AdvectionDiffusion name AdvectionDiffusion class $EXPN;
fem update field $nj_aA substitute constant 1; #ratio of a:A
fem update field 6 substitute constant 0;


#=================================
## FIRST BREATH
# Inspiration
$T_start = 0;
$T_end = $T_insp;
fem solve to $T_end delta_t $dt class $INSP;
fem export node;end_insp as end_insp field;
fem export elem;end_insp as end_insp field;
fem update solution class $EXPN substitute solution class $INSP;

# Expiration
$T_start = $T_end;
$T_end = $T_end + $T_expn;
fem update motion inlet_flow -$flow_in dv_from $INSP class $EXPN;
fem solve restart from $T_start to $T_end delta_t $dt class $EXPN;
fem update solution class $INSP substitute solution class $EXPN;
fem export node;end_exp as end_exp field;
fem export elem;end_exp as end_exp field;

fem export node;solution as airways;
fem export node;terminal as terminal nodes TERMINAL;
fem export elem;solution as airways;



# ## SUBSEQUENT BREATHS
# $NUM_BREATHS = 5;
# for $breath (2..$NUM_BREATHS){
#    # Inspiration
#    $T_start = $T_end;
#    $T_end = $T_end + $T_insp;
#    fem update solution class $INSP substitute solution class $EXPN;
#    fem update solution node 1 class $INSP substitute constant $conc;
#    fem update motion inlet_flow $flow_in dv_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;$outdir.$outfile."_end_insp" as end_insp field;
#    }
#    # Expiration
#    $T_start = $T_end;
#    $T_end = $T_end + $T_expn;
#    fem update solution class $EXPN substitute solution class $INSP;
#    fem update motion inlet_flow -$flow_in dv_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;$outdir.$outfile."_end_exp" as end_exp field;
#    }
# }



Additional testing commands:

fem list nodes

Files used by this example are:

Name                       Modified     Size

example_9a1.com 07-Apr-2011 4.3k  18-Feb-2011 95 AdvectionDiffusion.ipbase 18-Feb-2011 893 AdvectionDiffusion.ipequa 18-Feb-2011 1.8k AdvectionDiffusion.ippara 18-Feb-2011 5.8k AdvectionDiffusion.ipsolv 18-Feb-2011 2.5k Draw.com 18-Feb-2011 1.2k N2O2.ipmate 18-Feb-2011 288 cmiss_test.out 18-Feb-2011 9.9k draw.com 07-Apr-2011 1.4k radius.ipfiel 18-Feb-2011 5.4k test_geom.ipelem 18-Feb-2011 6.0k test_geom.ipnode 18-Feb-2011 7.5k test_output.com 18-Feb-2011 15 ventilation.ipfiel 18-Feb-2011 24k

Download the entire example:

Name                      Modified     Size

examples_9_9a_9a1.tar.gz 08-Apr-2011 30k

Testing status by version:

StatusTestedReal time (s)
i686-linux
cmFailureSun Mar 6 00:05:14 20160
last breakSun Feb 20 01:46:00 20111
cm-debugFailureSun Mar 6 00:05:15 20161
last breakSat Feb 19 01:36:00 20111
x86_64-linux
cmFailureSun Mar 6 00:01:25 20160
last breakSun Feb 20 01:10:00 20112
cm-debugFailureSun Mar 6 00:01:25 20160
last breakSat Feb 19 01:08:00 20110

Testing status by file:


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

Input last modified: Thu Apr 7 18:44:40 2011


CMISS Help / Examples / 9 / 9a / 9a1