sandbox/M1EMN/Exemples/viscolsqrt_noSV.c

    “full Huppert problem” 1D

    Equations

    We look at the following problem the the viscous collapse with hydrostatic pressure gradient of a lubricated film on a inclined plate of contant angle. This is as well a diffusive wave problem

    \displaystyle \frac{\partial}{\partial t} h + \frac{\partial}{\partial x} Q = 0,\text{ with } Q= \frac{h^3}3 (1 - \frac{\partial h}{\partial x})

    “Huppert’s problems”

    This is a mix of (Huppert first problem on a horizontal plate and second problem on an inclined plate:

    Huppert’s first problem is \displaystyle \frac{\partial}{\partial t} h = \frac{\partial}{\partial x}(\frac{h^3}{3} \frac{\partial}{\partial x} h ) whereas Huppert’s second one is \displaystyle \frac{\partial}{\partial t} h + h^2 \frac{\partial}{\partial x} h =0

    so that we have both effects here. This is refered as “Spreading of a shallow mass on an incline” by Mei. See bibliography

    Numerics

    We do not use here the shallow water solver (as in the two links below). We do a splitting to deal with the two parts of Q.

    set size ratio .3
set samples 9 
set label "h i-1" at 1.5,3.1
set label "h i" at 2.5,3.15
set label "h i+1" at 3.5,2.5
set xtics ("i-2" 0.5, "i-1" 1.5, "i" 2.5,"i+1" 3.5,"i+2" 4.5,"i+3" 5.5)
set arrow from 2,1 to 2.5,1
set arrow from 3,1 to 3.5,1
set label "Q i" at 2.1,1.25
set label "Q i+1" at 3.1,1.25

set label "x i-1/2" at 1.5,0.25
set label "x i" at 2.4,0.25
set label "x i+1/2" at 3.,0.25

set label "x"  at 0.5,2+sin(0) 
set label "x"  at 1.5,2+sin(1)
set label "x"  at 2.5,2+sin(2) 
set label "x"  at 3.5,2+sin(3) 
set label "x"  at 4.5,2+sin(4)  
f(x) = (2+sin(x))*(x<=4)
p[-1:7][0:4] f(x) w steps not,f(x) w impulse not linec 1
    (script)

    (script)

    The full problem is solved with a splitting Q=Q_c+ Q_d, each part has a different behavior. \displaystyle Q_c=\frac{h^3}3 \;\text{ and }\; Q_d= -\frac{h^3}3 \frac{\partial h}{\partial x} where subscript “c” refers to the “advective” part with a celerity c_0=h^2:

    \displaystyle \frac{\partial}{\partial t} h + c_0 \frac{\partial}{\partial x} h = 0 \;\;\text{ with }\;\; c_0 = \partial Q_c/\partial h

    and subscript “d” refers to the “diffusive” part of non linear viscosity \frac{h^3}{3}: \displaystyle \frac{\partial}{\partial t} h = - \frac{\partial}{\partial x} Q_d

    The first leads to shocks, the second is diffusive

    Code

    Mandatory declarations, libraries and tables. Note we use the face vectors.

    #include "grid/cartesian1D.h"
#include "run.h"

scalar h[];
scalar zb[];

face vector  hx[];
face vector Q[];
face vector c0[];
face vector Qc[];
face vector Qd[];

double alpha;
double X0,L0;
double DT;
double DD;

    Main with definition of parameters :

    Note that we run the code twice, the first run is without diffusion DD=0

    int main() 
 {
  L0 = 15;
  X0 = -1;
  N = 200;
  DT = (L0/N)*(L0/N)/20;
  DD = 0;
  run();
  DD = 1 ;
  run();
}

    Boundary conditions

    h[left] =  neumann (0);
hx.n[left] = neumann (0);
Q.x[left] = dirichlet(0);

    Initial elevation: a " quare" of surface 2

    event init (t = 0) 
 {
  foreach()
    h[] = (fabs(x-1)<1);
    
  foreach_face()
   {   
    Q.x[]=0;
    c0.x[]=1;
  }
}

    print data in stdout

    event printdata ( t+=5; t <= 100)
 {
  foreach()
    fprintf (stdout, "%g %g %g %g %g \n", x, h[], Q.x[], t, DD);
  fprintf (stdout, "\n");
}

    Integrating the domaine over time, updating the elevation at each time step.

    event integration (i++)
 {
    double dt = DT;

    finding the good next time step

        dt = dtnext (dt);

    the flux Q(h) : split in two fluxes.

    1- The avective part

    First Q_c the advective part gives a estimation of the height say h^*,

    \displaystyle \frac{h^* - h^n}{\Delta t} = -\frac{1}{\Delta x} (Q^n_{ai} - Q^n_{ai-1}) with Q_a obtained from Q_c with a correction

    \displaystyle Q_{ai} = \frac{Q_{ci} + Q_{ci-1}}{2} - c \frac{h_j - h_{j-1}}{2}\text{ with } c_0 = \frac{\partial Q_{c}}{\partial h}

        foreach_face()
      {
      double hm = ((h[0,0] + h[-1,0])/2);
      Q.x[] =  (pow(hm,3)/3) ;
      c0.x[] = pow(hm,2) *(1-DD);
      Qc.x[] = (Q.x[]+Q.x[-1])/2.  - (c0.x[] *(h[]-h[-1])/2);
      }
  
    foreach()
        h[] = h[] - dt*( Qc.x[1,0] - Qc.x[0,0] )/Delta;

    2- The diffusion part : calculate the next step elevation h^{n+1}, n: the temporal index.

    \displaystyle \frac{h^{n+1} - h^*}{\Delta t} = -\frac{1}{\Delta x} (Q^*_{di} - Q^*_{di-1})

        foreach_face()
     {
      hx.x[] = ((h[0,0] - h[-1,0] )/Delta);
      double hm = ((h[0,0] + h[-1,0])/2);
      Qd.x[] = -(pow(hm,3)/3) *  hx.x[];
    }

    final upadte, note the DD switch.
    the first run is without diffusion DD=0, the second is with: DD=1

        foreach()
      h[] = h[] - dt*DD*( Qd.x[1,0] - Qd.x[0,0] )/Delta;
    
}

    Plots

    Solution of the full problem (this is as well a diffusive wave problem).

    reset
set xlabel "x"
set size ratio .25
p'out' u 1:($5==1?$2:NaN)w l
    (script)

    (script)

    Solution of only the advected problem (second problem of Huppert, inclined plate, this is as well a kinetic wave problem); swith DD=0.

    reset
set xlabel "x"
set size ratio .25
p'out' u 1:($5==0?$2:NaN)w l
    (script)

    (script)

    Both superposed to see differences

    reset
set xlabel "x"
set size ratio .25
p'out' u 1:($5==0?$2:NaN) t 'adv' w l,'' u 1:($5==1?$2:NaN)t 'full' w l
    (script)

    (script)

    Plots of self similar solution

    Wave solution of form : \sqrt{x/t}=t^{-1/3} \sqrt{x/t^{1/3}}

    the solution is self similar at long times,

    reset
set xlabel "x"
set size ratio .25
p'out' u 1:($5==0?$2:NaN)w l,'' u 1:($5==1?$2:NaN)w l,'' u 1:(sqrt($1/$4)) w l
    (script)

    (script)

    the solution is self similar at long times, we plot in the selfsimilar space, we shift solution of the full problem (diffusive wave) to see better the differences induced by the addition of pressure gradient in the kinetic wave

    reset
set xlabel "x"
set size ratio .25
p[-1:4]'out' u ($1/$4**.33333333+$5):($2*$4**.3333333) not w l,sqrt(x),sqrt(x-1)
    (script)

    (script)

    Links

    Bibliography