#Example_952 Simulating water and heat transfer in a lung model
# 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;
}
#Parameters for mesh geometry
$RADIUS = 1;
# Names of major bronchi, for element groups
@bronchus = ("trachea","right_main_bronchus","left_main_bronchus","right_lower_bronchus","right_upper_lobe_bronchus","left_upper_lobe_bronchus","left_lower_lobe_bronchus","right_middle_lobe_bronchus","right_lower_lobe_bronchus");
# Start elements, end element,radius, and volume below for major bronchi
@data = ([1,128,8.0,0],[129,192,5.55,0],[193,256,6.0,0],[257,320,4.45,0],[321,384,3.65,1.5e4],[449,512,3.75,7.5e3],[385,448,4.0,2.0e4],[513,544,2.6,1.5e4],[545,576,3.2,2.0e4]);
$start = 0;
$end = 1;
$size = 2;
$volu = 3;
#Parameters for problem solution
$INSPIRATION = 1; #class number for inspiration
$EXPIRATION = 2; #class number for expiration
$dt = 0.05; #size of time step (s)
$T_insp = 0.2; #duration of inspiration (s)
$T_expn = 0.2; #duration of expiration (s)
$NUM_EXPORTS=1;
$time = 0.0; #initial time (s)
$flow_in = 0.1; #volume of inspiration (litres)
$flow_out = -0.1; #volume of expiration (litres)
$FIXED_NODE = 1; #node to apply fixed boundary condition (inspiration)
$FLUX_NODE = 1; #node to apply flux boundary condition (expiration)
$NUM_BREATHS = 2; #number of breaths (inspiration + expiration)
$N=0;
fem de param;r;$example/fullung;
fem de coor;r;$example/fullung;
fem de base;r;$example/fullung; #linear basis function, 4 gp
fem de node;r;$example/fullung; #define model nodes
fem de elem;r;$example/fullung 1d; #define model elements
# Set up default airway mesh structure. i.e. create
# fields with versions at the bifurcation points
fem de mesh;c airway as conducting_airways field;
# Evaluate the generations, Horsfield orders, and Strahler orders
fem evaluate ordering;
# Assign radii to the element group, using Weibel data for generation
fem update mesh geometry elem conducting_airways radius_field $RADIUS anatomical weibel alveolar_field 2;
$i = 0;
foreach $bronchus(@bronchus){
# Group the elements in each major bronchus
fem group elem "$data[$i][$start]..$data[$i][$end]" as $bronchus;
# Assign radii for the major bronchi using Horsfield values
fem up mesh geom elem $bronchus radius_field $RADIUS const $data[$i][$size];
}
# Add lumped parameter models to each peripheral airway
fem de;add mesh;c airway lumped_parameter as respiratory field radius_field 1 alveolar_field 2;
for my $i (4..8){
fem group elem all parent $data[$i][$end] as dummy_a;
# Group terminal branches in the lobe == lumped parameter/respiratory elements
fem group elem dummy_a384 by terminal as dummy_b;
# Give the lumped parameter models volumes specific to the lobe
fem up mesh geom elem dummy_b radius $RADIUS lumped volume $data[$i][$volu];
}
#DEFINE THE PROBLEM SETUP FOR INSPIRATION
fem de equa;r;$example/fullung class $INSPIRATION;
fem de mate;r;$example/fullung class $INSPIRATION;
fem de moti;c lung elem respiratory uniform flow_field 3 inlet_flow $flow_in class $INSPIRATION;
fem de init;c lung fixed node $FIXED_NODE temperature 31.3 rh 100 core 37 class $INSPIRATION;
fem de solve;r;$example/fullung class $INSPIRATION;
#DEFINE THE PROBLEM SETUP FOR EXPIRATION
fem de equa;r;$example/fullung class $EXPIRATION;
fem de mate;r;$example/fullung class $EXPIRATION;
fem de moti;c lung elem respiratory uniform flow_field 3 inlet_flow $flow_out class $EXPIRATION;
fem de init;c lung flux node $FLUX_NODE class $EXPIRATION;
fem de solve;r;$example/fullung class $EXPIRATION;
$restart = "";
$time_end = 0.0;
for (my $j = 0; $j < $NUM_BREATHS; $j++){
for (my $i = 0; $i < $NUM_EXPORTS; $i++){
$time_end = $time_end + $T_insp/$NUM_EXPORTS;
print "FEM solve $restart to $time_end delta_t $dt class $INSPIRATION\n";
fem solve $restart to $time_end delta_t $dt class $INSPIRATION;
$N++;
$restart = "restart";
}
fem update solution class $EXPIRATION substitute solution class $INSPIRATION;
for (my $i = 0; $i < $NUM_EXPORTS; $i++){
$time_end = $time_end + $T_expn/$NUM_EXPORTS;
print "FEM solve $restart to $time_end delta_t $dt class $EXPIRATION\n";
fem solve $restart to $time_end delta_t $dt class $EXPIRATION;
$N++;
}
fem update solution class $INSPIRATION substitute solution class $EXPIRATION;
fem update solution class $INSPIRATION nodes $FIXED_NODE substitute constant (31.3+273.15);
}