Comment Constructs slowly-varying lubrication model of thin
film 3D fluid flows on arbitrarily curved 2D substrates.
Allows for large changes in film thickness.  Incorporates
surface tension, gravity, inertia, slip.  The coordinate 
system is (x,z,y), velocity (u,w,v), scale factors 
(h1,h2,1).

Last updated 16 Sept 2003;

% \section{Computer algebra derives the model}
% \label{3dp}
% 
% The \textsc{reduce}\footnote{At the time of writing, information about
% \textsc{reduce} was available from Anthony C.~Hearn, RAND, Santa
% Monica, CA~90407-2138, USA. \protect\url{mailto:reduce@rand.org}}
% program which performs the derivation of centre manifold model is
% listed below.  A little explanation is usefully given first.  Observe
% that program variables are typeset in \verb|teletype| font.
% \begin{itemize}
%     \item The physical coordinate system within the program is 
%     $(x_1,x_2,y)=(\verb|x|,\verb|z|,\verb|y|)$ with scale factors 
%     $(h_1,h_2,h_3)=(\verb|h1|,\verb|h2|,1)$\,, velocity field 
%     $(u_1,u_2,v)=(\verb|u|,\verb|w|,\verb|v|)$ and pressure 
%     $p=\verb|p|$\,.  Scale factors of the substrate are $m_1=\verb|m1|$ and 
%     $m_2=\verb|m2|$\,, whereas those evaluated on the free surface are 
%     $\tilde h_1=\verb|hh1|$ and~$\tilde h_2=\verb|hh2|$\,.  
%     
%     Substrates with specific geometries, as discussed in 
%     \S\ref{Sspec}, may be derived simply by coding information 
%     about their curvatures~$k_i$ and scale factors~$m_i$ as shown for 
%     the three cases discussed in \S\ref{Sspec}.
% 
%     \item Expressions are written in terms of a stretched coordinate 
%     system~$X=x_1$, $Z=x_2$, $Y=y/\eta$, $T=t$ so that the free 
%     surface of the fluid film is simply~$Y=1$\,.  In the program we use 
%     $(X,Z,Y,T)=(\verb|xs|,\verb|zs|,\verb|ys|,\verb|ts|)$.  Then 
%     spatio-temporal derivatives transform by the chain rule
%     \begin{displaymath}
%         \D{x_1}{}=\D X{}-\frac{Y\eta_X}{\eta}\D Y{}\,,\quad
% 	\D{x_2}{}=\D Z{}-\frac{Y\eta_Z}{\eta}\D Y{}\,,
%     \end{displaymath}
%     \begin{displaymath}
%         \D{t}{}=\D T{}-\frac{Y\eta_T}{\eta}\D Y{}\,,\quad
% 	\mbox{and}\quad
% 	\D{y}{}=\frac{1}{\eta}\D Y{}\,,
%     \end{displaymath}
%     as coded.
% 
%     \item The amplitudes of the model are $(\eta,\bu _1,\bar 
%     u_2)=(\verb|h|,\verb|uu|,\verb|ww|)$ with \verb|h(m,n)| denoting 
%     the spatial derivative $\frac{\partial^{m+n} \eta}{\partial X^m 
%     \partial Z^n}$ and similarly for \verb|uu(m,n)| and 
%     \verb|ww(m,n)|.  The evolution of these quantities is given by
%         \begin{displaymath}
%         \D t\eta=\verb|gh|\,,\quad
% 	\D t{\bu _1}=\verb|gu|\,,\quad\mbox{and}\quad
% 	\D t{\bu _2}=\verb|gw|\,.
%     \end{displaymath}
% \end{itemize}
% The correctness of the results of this program depend only upon the 
% correct coding of the physical equations.  The algebraic machinations 
% are repeated until the residual of the fluid differential equations 
% and boundary conditions are zero to the requisite order.


% improve printing appearance
on div; off allfac; on revpri;
factor d,h,uu,ww,re,gravity,betag,gr,we,k1,k2,sk;
% use operator h(m,n) to denote df(h,x,m,z,n)
operator h;
depend h,xs,zs,ts;
let { df(h(~m,~n),xs) => h(m+1,n)
     ,df(h(~m,~n),zs) => h(m,n+1)
     ,df(h(~m,~n),xs,zs) => h(m+1,n+1)
     ,df(h(~m,~n),ts) => df(gh,xs,m,zs,n) };
% principal curvatures of the substrate
operator k1;
depend k1,xs,zs;
let { df(k1(~m,~n),xs) => k1(m+1,n)
     ,df(k1(~m,~n),zs) => k1(m,n+1) };
operator k2;
depend k2,xs,zs;
let { df(k2(~m,~n),xs) => k2(m+1,n)
     ,df(k2(~m,~n),zs) => k2(m,n+1) };
operator uu; operator ww;
depend uu,xs,zs,ts;
depend ww,xs,zs,ts;
let { df(uu(~m,~n),xs) => uu(m+1,n)
     ,df(ww(~m,~n),xs) => ww(m+1,n)
     ,df(uu(~m,~n),zs) => uu(m,n+1)
     ,df(ww(~m,~n),zs) => ww(m,n+1)
     ,df(uu(~m,~n),ts) => df(gu,xs,m,zs,n)
     ,df(ww(~m,~n),ts) => df(gw,xs,m,zs,n) };
% use stretched coords: ys=y/h(x,z,t), xs=x, zs=z, ts=t
depend xs,x,y,z,t;
depend zs,x,y,z,t;
depend ys,x,y,z,t;
depend ts,x,y,z,t;
let { df(~a,x) => df(a,xs)-ys*h(1,0)/h(0,0)*df(a,ys)
    , df(~a,z) => df(a,zs)-ys*h(0,1)/h(0,0)*df(a,ys)
    , df(~a,t) => df(a,ts)-ys*gh/h(0,0)*df(a,ys)
    , df(~a,y) => df(a,ys)/h(0,0)
    };
% some abbreviations with appropriate scalings
kx:=d*k1(0,0)$ kz:=d*k2(0,0)$
eta:=h(0,0)$ etax:=d*h(1,0)$ etaz:=d*h(0,1)$
% scale factors of coord system: substrate, fluid & surface
depend m1,xs,zs;
depend m2,xs,zs;
h1:=m1*(1-kx*eta*ys)$
h2:=m2*(1-kz*eta*ys)$
hh1:=sub(ys=1,h1)$
hh2:=sub(ys=1,h2)$
% gravitational body force
depend g1,xs,zs;
depend g2,xs,zs;
depend gy,xs,zs;
let { df(g1,xs) => -rh2*df(h1,z)*g2+m1*k1(0,0)*gy
    , df(g1,zs) => rh1*df(h2,x)*g2
    , df(g2,zs) => -rh1*df(h2,x)*g1+m2*k2(0,0)*gy
    , df(g2,xs) => rh2*df(h1,z)*g1
    , df(gy,xs) => -m1*k1(0,0)*g1
    , df(gy,zs) => -m2*k2(0,0)*g2
    };
% computes mean across fluid layer
operator mean; linear mean;
let{ mean(ys^~n,ys)=> 1/(n+1)
    ,mean(ys,ys)=>1/2
    ,mean(1,ys)=>1 };
% solves -df(p,ys)=rhs s.t. sub(ys=1,p)=0
operator psolv; linear psolv;
let {psolv(ys^~n,ys) => (1-ys^(n+1))/(n+1)
    ,psolv(ys,ys) => (1-ys^2)/2
    ,psolv(1,ys) => (1-ys) };
% solves df(u,ys,2)=rhs s.t. sub(ys=0,u)=0 & mean(u)=0
% and df(v,ys,2)=rhs s.t. sub(ys=0,u)=0 & mean(v)=0
operator usolv; linear usolv;
let {usolv(ys^~n,ys) => (ys^(n+2)-2*ys/(n+3))/(n+2)/(n+1)
    ,usolv(ys,ys) => (ys^3-2*ys/4)/3/2
    ,usolv(1,ys) => (ys^2-2*ys/3)/2 };

% may derive model in special geometries if activated
if 0 then % flat substrate in Cartesian
let { k1(~p,~q)=>0, m1=>1
    , k2(~p,~q)=>0, m2=>1 } 
else if 0 then % cylindrical substrate of radius rad
% x1 axial distance & x2 is angle around
% sk=+1 is flow outside cyl, sk=-1 is flow inside cyl
let { k1(~p,~q)=>0 , m1=>1
    , k2(0,0)=>-sk/rad , sk^2=>1
    , k2(~p,~q)=>0 when p+q>0 , m2=>rad } 
else if 1 then begin % substrate depth(zs) below reference
% s=x1 longitudinal distance & r=x2 is cross stream location
% dd=df(depth,x2), m2=\sqrt{1+dd^2}, k2=-dd'/m2^3
depend dd,zs; 
let { k1(~p,~q)=>0, m1=>1
    , k2(~p,~q)=>0 when p>0
    , df(m2,xs)=>0, df(m2,zs)=>-dd*m2^2*k2(0,0)
    , df(dd,zs)=>-m2^3*k2(0,0) 
    , dd^2=>m2^2-1 } end
else if 0 then % flat substrate in polar coords
% x1 radial distance & x2 is angle around
let { k1(~p,~q)=>0 , m1=>1
    , k2(~p,~q)=>0 , m2=>xs } 
else if 0 then % outside spherical coords x1=colat x2=long
let { k1(0,0)=>-1/rad 
    , k1(~p,~q)=>0 when p+q>0 , m1=>rad
    , k2(0,0)=>-1/rad 
    , k2(~p,~q)=>0 when p+q>0 , m2=>rad*sin(xs) } ;

% centre subspace and null evolution
u:=d*2*ys*uu(0,0); w:=d*2*ys*ww(0,0); v:=0; p:=0;
gh:=0; gu:=0; gw:=0;
% and linear approx to various reciprocals
% rh1, rh2 are the approx of 1/h1, 1/h2
rh1:=1/m1$ rh2:=1/m2$ rhh1:=rh1$ rhh2:=rh2$

% Use d to count number of uu, ww and derivatives of x & z
% Throw away this order or higher in uu, ww, d/dx & d/dz
let {d^4=>0, gam^8=>0};  
gravity:=d^2*gr; % gravity is the coefficient/order of gravity terms
% Iterate until converges
it:=0$
repeat begin
   write it:=it+1;
   % approximate mean curvature of the surface
   curv:=sub(ys=1,kx+kz +d*rh1*rh2*(df(h2*rh1*etax,x)
             +df(h1*rh2*etaz,z)) +(kx^2+kz^2)*eta);
   % free-surface deviatoric stress (Batchelor, p600)
   txx:=sub(ys=1,2*d*(rhh1*df(u,x)+rhh1*rhh2*df(h1,z)*w
                 +rhh1*v*df(h1,y)));
   txz:=sub(ys=1,d*(rhh2*df(u,z)+rhh1*df(w,x)
                 -rhh1*rhh2*(df(h1,z)*u+df(h2,x)*w)));
   txy:=sub(ys=1,hh1*df(rh1*u,y)+rhh1*d*df(v,x));
   tzz:=sub(ys=1,2*d*(rhh2*df(w,z)+rhh1*rhh2*df(h2,x)*u
                 +rhh2*v*df(h2,y)));
   tzy:=sub(ys=1,hh2*df(rh2*w,y)+rhh2*d*df(v,z));
   tyy:=sub(ys=1,2*df(v,y));
   % omega=curl(q) => del^2q=-curl(omega) (Batchelor p599)
   o1:=rh2*(d*df(v,z)-df(h2*w,y));
   o2:=rh1*(df(h1*u,y)-d*df(v,x));
   oy:=rh1*rh2*(d*df(h2*w,x)-d*df(h1*u,z));
   % vertical momentum & normal surface stress
   begin scalar veq,ns;
   veq:=re*( df(v,t)+d*u*rh1*df(v,x)+v*df(v,y)
            +d*w*rh2*df(v,z)+d*m1*k1(0,0)*rh1*u^2
        +d*m2*k2(0,0)*rh2*w^2 )
        +df(p,y)-gravity*gy+rh1*rh2*d*(df(h2*o2,x)-df(h1*o1,z));
   ns:=(hh2^2*etax^2*txx+2*hh1*hh2*etax*etaz*txz
       -2*hh1*hh2^2*etax*txy -2*hh1^2*hh2*etaz*tzy
       +hh1^2*etaz^2*tzz+hh1^2*hh2^2*tyy)
       -sub(ys=1,p)*((hh2*etax)^2+(hh1*etaz)^2+(hh1*hh2)^2)
       -we*curv*((hh2*etax)^2+(hh1*etaz)^2+(hh1*hh2)^2);
   ok:=if (veq=0)and(ns=0) then 1 else 0;
   pd:=eta*psolv(veq,ys)+ns*(rhh1*rhh2)^2;
   p:=p+pd;
   end;
   % x1-momentum & tangential stress on FS
   begin scalar ueq,ts1,gud,ubed;
   ueq:=re*( df(u,t)+d*u*rh1*df(u,x)+v*df(u,y)
            +d*w*rh2*df(u,z)+d*w*rh1*rh2*(u*df(h1,z)
        -w*df(h2,x))-d*m1*k1(0,0)*rh1*v*u )
        +rh1*d*df(p,x)-gravity*g1+rh2*(d*df(oy,z)-df(h2*o2,y));
   ts1:=((hh1^2*hh2-hh2*etax^2)*txy
         +hh1*hh2*etax*(tyy-txx)-hh1^2*etaz*txz
         -hh1*etax*etaz*tzy )*rhh1^2*rhh2
         -(1-gam)*sub(ys=1,u)/eta;
   ubed:=sub(ys=0,-u);
   ok:=if ok and(ueq=0)and(ts1=0)and(ubed=0) then 1 else 0;
   gud:=(-mean(ueq*ys,ys)-ts1/eta)*3/2;
   gu:=gu+gud/re/d;
   u:=u+usolv(ueq+2*ys*gud,ys)*eta^2+ubed;
   end;
   % x2-momentum & tangential stress on FS
   begin scalar weq,ts2,gwd,wbed;
   weq:=re*( df(w,t)+d*u*rh1*df(w,x)+v*df(w,y)
            +d*w*rh2*df(w,z)+d*u*rh1*rh2*(w*df(h2,x)
        -u*df(h1,z))-d*m2*k2(0,0)*rh2*v*w )
        +rh2*d*df(p,z)-gravity*g2+rh1*(df(h1*o1,y)-d*df(oy,x));
   ts2:=((hh1*hh2^2-hh1*etaz^2)*tzy
       +hh1*hh2*etaz*(tyy-tzz)-hh2*etax*etaz*txy
       -hh2^2*etax*txz)*rhh1*rhh2^2
       -(1-gam)*sub(ys=1,w)/eta;
   wbed:=sub(ys=0,-w);
   ok:=if ok and(weq=0)and(ts2=0)and(wbed=0) then 1 else 0;
   gwd:=(-mean(weq*ys,ys)-ts2/eta)*3/2;
   gw:=gw+gwd/re/d;
   w:=w+usolv(weq+2*ys*gwd,ys)*eta^2+wbed;
   end;
   % continuity & bed
   begin scalar ceq;
   ceq:=-( d*df(h2*u,x)+d*df(h1*w,z)+df(h1*h2*v,y))/m1/m2;
   ok:=if ok and(ceq=0) then 1 else 0;
   v:=v+eta*int(ceq,ys);
   end;
   % uu and ww represent weighted mean velocity
   begin scalar ueq,weq;
   ueq:=uu(0,0)*d-mean(u*h2,ys)/m2;
   weq:=ww(0,0)*d-mean(w*h1,ys)/m1;
   ok:=if ok and(weq=0)and(ueq=0) then 1 else 0;
   u:=u+2*ys*ueq;
   w:=w+2*ys*weq;
   end;
   % refine approx of h1 and h2
   h1eq:=rh1*h1-1$ rh1:=rh1-h1eq/m1; rhh1:=sub(ys=1,rh1);
   h2eq:=rh2*h2-1$ rh2:=rh2-h2eq/m2; rhh2:=sub(ys=1,rh2);
   write ok:=if ok and(h1eq=0)and(h2eq=0) then 1 else 0;
   % kinematic BC on FS
   gh:=sub(ys=1,v-u*rh1*etax-w*rh2*etaz);
   showtime;
end until ok or(it>19);

if 0 then begin % get shear stress on the substrate
    txy:=sub(ys=0,hh1*df(rh1*u,y)+rhh1*d*df(v,x))$
    tzy:=sub(ys=0,hh2*df(rh2*w,y)+rhh2*d*df(v,z))$
    gam:=1;
    on rounded;
    print_precision 4;
    coeffn(txy,d,1);
    coeffn(tzy,d,1);
    coeffn(txy,d,2);
    coeffn(tzy,d,2);
    coeff(gu,d);
    coeff(gw,d);
end;

end;