Advection.proj
This is a pure advection or wave equation in a tube. What comes in at the left
hand boundary exits at the right hand boundary.
There are two independent variables, t in seconds, where t stands
for time, and x in centimeters, where x is a spatial variable.
The tube is parameterized by length L, volume V, and flow F.
The velocity in the tube, U, is given by U=F*L/V. The
transit time, Tmean (time for the input concentration to appear at the
outflow is given by Tmean = V/F.
The left hand Dirichlet boundary condition is given as
-U*C = -U*Cin.
The right hand Neumann boundary condition, C:x=0, is, strictly speaking,
extraneous, insofaras there is only one derivative in the x space. However,
JSim requires two boundary conditions for PDE's.
The result is the advection of the wave form (Lagged-normal curve
representing a bolus injection) in the downstream direction.
(see Figures 1 and 2).
import nsrunit; unit conversion on;
// SHORT DESCRIPTION
// Model for the flow of a substance in a tube.
math Advection {
realDomain t sec; t.min=0; t.max=15; t.delta=0.1;
real L = 0.1 cm, // Length of tube
Ngridx = 51; // Number of grid points spatially
realDomain x cm; x.min=0; x.max=L; x.ct=Ngridx;
real C(x,t) mM; // Concentration of substance
extern real Cin(t) mM; // Inflow Concentration
extern real C0(x) mM; // Initial Distribution of C
real Cout(t) mM, // OutflowConcentration
F = 1 ml/min, // Flow rate
V = 0.05 ml, // Volume of tube
U cm/sec, // Flow Velocity
Tmean sec; // Time for inflow to appear as outflow
U = F*L/V;
Tmean = V/F;
// Initial Condition
when(t=t.min) {C(x,t)=if(x=x.min) Cin else C0;}
// Boundary Conditions
when (x=x.min) { -U*C = -U*Cin;} // Left Hand BC
when (x=x.max) { C:x = 0; Cout = C;} // Right Hand BC
// Partial Differential Equation
C:t = -U*C:x ;
}