COMMENT find the flow in an arbitrarily curving pipe
Simultaneously determine the shear dispersion in such a flow:
eps=magnitude of curvature terms,
del=magnitude of axial derivatives and of torsion.
;
on div; off allfac; on revpri;
factor eps,del,cos,sin,df,kap,tau;
%
depend kap,s; % curvature
depend tau,s; % torsion
% local coordinate system is (s,r,th), tp=th+phi(s)
depend tp,s,th;
let { df(tp,s)=>-tau, df(tp,th)=>1 };
hs:=1-eps*kap*r*cos(tp);
hr:=1;
ht:=r;
% trigonometry rules OK
let { sin(~a)*cos(~b) => (sin(a+b)+sin(a-b))/2
    , cos(~a)*cos(~b) => (cos(a-b)+cos(a+b))/2
    , sin(~a)*sin(~b) => (cos(a-b)-cos(a+b))/2
    , cos(~a)^2       => (1+cos(2*a))/2
    , sin(~a)^2       => (1-cos(2*a))/2
    };
% get the streamwise velocity (mult rhs by r^2 to apply)
depend r,rt;depend tp,rt;
operator usol; linear usol;
let { usol(r^~m*cos(~n*tp),rt) => (r^m-r^n)*cos(n*tp)/(m^2-n^2)
    , usol(r^~m*sin(~n*tp),rt) => (r^m-r^n)*sin(n*tp)/(m^2-n^2)
    , usol(r^~m*cos(tp),rt) => (r^m-r)*cos(tp)/(m^2-1)
    , usol(r^~m*sin(tp),rt) => (r^m-r)*sin(tp)/(m^2-1)
    , usol(r^~m,rt) => (r^m-1)/m^2
    };
% mean over a cross section (mult by r to use)
operator mean; linear mean;
let { mean(r^~m*cos(~n),rt) => 0
    , mean(r^~m*sin(~n),rt) => 0
    , mean(r^~m,rt) => 2/(m+1)
    , mean(r,rt) => 1
    };
%
% get the radial velocity (mult by r^2 to use)
operator vsol; linear vsol;
let { vsol(r^~m*cos(~n*tp),rt) => 
      (r^m-((1+n-m)*r^(n-1)+(1-n+m)*r^(n+1))/2)*cos(n*tp)
      /((m-1)^2-n^2)/((m+1)^2-n^2)
    , vsol(r^~m*sin(~n*tp),rt) => 
      (r^m-((1+n-m)*r^(n-1)+(1-n+m)*r^(n+1))/2)*sin(n*tp)
      /((m-1)^2-n^2)/((m+1)^2-n^2)
    , vsol(r^~m*cos(tp),rt) => (r^m-((2-m)+m*r^2)/2)*cos(tp)
                               /((m-1)^2-1)/((m+1)^2-1)
    , vsol(r^~m*sin(tp),rt) => (r^m-((2-m)+m*r^2)/2)*sin(tp)
                               /((m-1)^2-1)/((m+1)^2-1)
    , vsol(r^~m,rt) => (r^m-r)/(m-1)^2/(m+1)^2
    };
%
% swirling velocity
operator wsol; linear wsol;
let { wsol(r^~m*cos(~n),rt) => 0
    , wsol(r^~m*sin(~n),rt) => 0
    , wsol(r^~m,rt) => (r^m-r)/(m^2-1)
    };
%
% radial pressure gradient
operator psol; linear psol;
let { psol(r^~m*cos(~n),rt) => 0
    , psol(r^~m*sin(~n),rt) => 0
    , psol(r*cos(~n),rt) => 0
    , psol(r*sin(~n),rt) => 0
    , psol(r^~m,rt) => r^m/m
    , psol(1,rt) => 1
    };
%
% to solve for the tracer concentration (mult by r^2 to use)
operator csol; linear csol;
let { csol(r^~m*cos(~n*tp),rt) => (r^m-m/n*r^n)*cos(n*tp)/(m^2-n^2)
    , csol(r^~m*sin(~n*tp),rt) => (r^m-m/n*r^n)*sin(n*tp)/(m^2-n^2)
    , csol(r^~m*cos(tp),rt) => (r^m-m*r)*cos(tp)/(m^2-1)
    , csol(r^~m*sin(tp),rt) => (r^m-m*r)*sin(tp)/(m^2-1)
    , csol(r^~m,rt) => (r^m-2/(m+2))/m^2
    };
%
% angular integration, but ditch secular terms !!!!
operator intt; linear intt;
let { intt(cos(~n*tp),tp) =>  sin(n*tp)/n
    , intt(sin(~n*tp),tp) => -cos(n*tp)/n
    , intt(cos(tp),tp) =>  sin(tp)
    , intt(sin(tp),tp) => -cos(tp)
    , intt(1,tp) => 0
    };

% straight tube limit of Poiseille flow, with mean velocity 1
u:=2*(1-r^2) +eps*3/2*kap*(r-r^3)*cos(tp);
v:=0;
w:=0;
% pressure = ps + p 
% = (local mean gradient) + (zero mean fluctuation)
ps:=-8; p:=0;
% concentration of tracer, mean c.
depend c,s,t;
let df(c,t)=>g;
cc:=c;
g:=0;
pe:=re*sc; % Peclet number = Reynolds * Schmidt
rh:=1; % approx reciprocal of axial scale factor hs

% let { df(kap,s)=>0, df(tau,s)=>0 }; % use for helical pipe
% re:=0; % use for Stokes flow
% sc:=81/10;  % use for heat in water at 15C
% sc:=71/100; % use for heat in air at 15C
% iterate until residuals are negligible
let { eps^3=>0, del^4=>0 }; % truncate the asymptotics
repeat begin 
    % reciprocal of scale factor
    eqr:=hs*rh-1;
    rh:=rh-eqr;
    % vorticity
    oms:=(df(r*w,r)-df(v,th))/r;
    omr:=(df(hs*u,th)-del*df(r*w,s))*rh/r;
    omt:=(del*df(v,s)-df(hs*u,r))*rh;
    % Navier-Stokes equation
    nss:=rh*(ps+del*df(p,s)) +(df(r*omt,r)-df(omr,th))/r
        +re*( del*u*df(u,s)*rh+v*df(u,r)+w*df(u,th)/r 
             +u*df(hs,r)*v*rh+u*df(hs,th)*w*rh/r );
    nsr:=df(p,r)    +(df(hs*oms,th)-del*df(r*omt,s))*rh/r
        +re*( del*u*df(v,s)*rh+v*df(v,r)+w*df(v,th)/r
             -w^2/r-df(hs,r)*u^2*rh );
    nst:=df(p,th)/r +(del*df(omr,s)-df(hs*oms,r))*rh
        +re*( del*u*df(w,s)*rh+v*df(w,r)+w*df(w,th)/r 
             -u^2*df(hs,th)*rh/r+v*w/r );
    % continuity equation
    cty:=(del*df(r*u,s)+df(r*hs*v,r)+df(hs*w,th))*rh/r;
    % find u and correction to mean pressure gradient
    ud:=usol(r^2*nss,rt);
    udm:=mean(r*ud,rt);
    u:=u+ud-udm*2*(1-r^2);
    ps:=ps+udm*8;
    % manipulate residuals
    eq1:=nsr-df(r^2*cty,r)/r^2;
    eq2:=df(nst,th)+df(df(r^2*cty,r)/r,r);
    eq3:=df(eq1,th,2)-df(r*eq2,r);
    % find v, w, p
    vd:= vsol(r^2*eq3,rt);
    wd:=-intt(r*cty+df(r*vd,r),tp);
    wd:=wd+wsol(r^2*(nst-df(df(r*wd,r)/r,r)+df(vd/r,th,r)),rt);
    pd:=-intt(r*(nst-df(df(r*wd,r)/r,r)+df(vd/r,th,r)),tp)
        -psol(r*nsr,rt);
    pdm:=mean(r*pd,rt);
    v:=v+vd;
    w:=w+wd;
    p:=p+pd-pdm; ps:=ps+pdm;
    % equation for tracer evolution
    ceq:= df(cc,t) +pe*( del*u*df(cc,s)*rh+v*df(cc,r)+w*df(cc,th)/r )
         -(del^2*df(rh*r*df(cc,s),s)+df(r*hs*df(cc,r),r)
            +df(hs/r*df(cc,th),th))*rh/r;
    cmean:=mean(r*hs*cc,rt)-c;
    % solvability & solution
    gd:=-mean(r*ceq,rt);
    g:=g+gd;
    cc:=cc+csol(r^2*(ceq+gd),rt)-cmean;
    showtime;
end until (eqr=0)and(nss=0)and(nsr=0)and(nst=0)
       and(cty=0)and(ceq=0)and(cmean=0);

% check subsiduary conditions
uwall:=sub(r=1,u);
vwall:=sub(r=1,v);
wwall:=sub(r=1,w);
umean:=mean(r*u,rt);
pmean:=mean(r*p,rt);
cwall:=sub(r=1,df(cc,r));
cflux:=mean(r*(u*pe*cc-del*rh*df(cc,s)),rt);
off nat;
out piped_o;
cflux:=cflux;
shut piped_o;

end;