sandbox/M1EMN/Exemples/darcysilo.c

    Solution of the potential flow though a small orifice in the bottom of a silo (by Matched Asymptotic Expansion)

    the problem :

    solution of the “Basic problem” in 2D, potential incompressible flow \displaystyle u =\frac{\partial p}{\partial x}, \;\; v = \frac{\partial p}{\partial y},\;\;\; \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}=0 to be solved in the upper half space 1>y>0 and -0.5 < x < 0.5 (the walls are in x =\pm 1/2 y>0, and \frac{\partial p}{\partial x}=0 on the walls) filled with a porous media, pressure is given at the top. At the bottom, there is a very small slit (- \varepsilon < x < \varepsilon, y=0) where pressure is constant: p=0. We take \varepsilon=.1. And there is no penetration for (|x|>.1, y=0), so normal velocity is \frac{\partial p}{\partial y}=0

    There is here a small parameter \varepsilon, we do a Matched Asymptotic Expansion using this small parameter.

    solution in complex plane for asymptotic approximations of the problem

    Inner Problem: Lamb

    First we look at small \varepsilon near the slit itself, so we rescale x=\varepsilon \tilde x, y=\varepsilon \tilde y, pressure stays at the same scale. The problem is now to solve the laplacian of pressure in the upper half plane \tilde y >0, with the slit (-c < \tilde {x} < c, (here c=1) \tilde y=0) where pressure is constant, say \tilde p=P_0. And there is no penetration on the wall for (| \tilde x | > c, \tilde y=0), so normal velocity is \frac{\partial \tilde p}{\partial \tilde y}=0

    Solution for the case of a slit with the upper half plane \tilde >0, full of fluid is proposed in Sneddon and in Lamb (c is the width of the slit, which is one half). It is \tilde x=c cosh(\tilde p) cos (\tilde \psi) and \tilde y=c sinh(\tilde p) sin (\tilde \psi)

    with \tilde z/c = cosh(F(\tilde z)) with \tilde z = \tilde x + i \tilde y and F=\tilde p + i \tilde psi a complex potential. iso pressure are ellipses (\tilde x/c/cosh(\tilde p))^2 + (\tilde y/c/sinh(\tilde p))^2=1. Far enough sinh(\tilde p) \simeq cosh(\tilde p) \simeq e^{\tilde p}/2, so that \displaystyle \tilde p = P_0 + log (\tilde x^2 + \tilde y^2) + log (2)

    Note that the P_0 is an extra constant that we add for convinioance and matching, the strict Lamb solution is with P_0=0 (this solution has been tested here darcyLambSneddon.c with P_0=0) far enough it is a source.

    Outer problem: Paterson

    We turn now to the external problem, for large x and y measured in \varepsilon.

    Hence we guess that we have to look at the problem of a source/sink Log(x+i y) between two walls in x=\pm 1.

    A king of method of images is used: we superpose an infinite number of sources/sinks (a =1 distance between one and the other). Solution of the problem (for a very small slit c \ll a) is \displaystyle \Sigma_{-\infty}^{\infty}Log(z - n a)

    The same problem exists for an infinite sum of vortices: (i \Sigma Log(z - n a)). It was solved by Karman (as cited by Péres 1936), and is in Paterson (1984). Thanks to Weierstrass Borel factorization: \displaystyle sin ( \pi z)/(\pi z) = \Pi_1^{\infty} (1 -z^2/n^2) we take this identity \displaystyle sin ( \pi z) = \pi z (1 -z^2) (1 -z^2/4 )... develop it \displaystyle sin (\pi z) = \pi z (1 -z)(1+z) (1 -z/2 ) (1+z/2) ... hence \displaystyle Log sin ( \pi z) = Log ( \pi ) + Log ( z) + Log ( 1 -z) + Log ( 1+z ) + Log ( 1+z/2 ) + Log ( 1- z/2 ) +... we change the order to recognize our sum \displaystyle Log sin ( \pi z) = Log ( \pi ) + Log ( z) + Log ( 1 -z) + Log ( 1+z ) + Log ( 2- z ) + Log(1/2) + Log ( 2+ z ) + Log(1/2) +... which is \displaystyle Log sin ( \pi z) = Log ( \pi ) + \Sigma_{n} Log(z - n ) + \Sigma_{n>0} 2 Log( 1/n)

    There is a “constant” which is fact infinite (Log(\pi)+2\Sigma_{n>0} Log( 1/n)), as Paterson says “except a constant which does not affect velocity”

    Final expression for complex potential \displaystyle Log(sin(\pi z))= \Sigma_{n=-\infty}^{n=\infty} Log(z - n )+( 2 \Sigma_{n=1}^{\infty} Log( 1/n)+ Log( \pi))

    Matching

    for small x and y close to the sink, the outer problem: \displaystyle Log(sin(\pi z))= Log \pi + Log z + ... this gives the outer pressure (a sink): \displaystyle p= Log \pi + Log r + ... but from Lamb inner problem, far from the slit, far away, cosh( \tilde p) \simeq sinh( \tilde p) \simeq e^{\tilde p}/2, hence iso pressure are \tilde p \simeq P_0+ log(2 \sqrt{\tilde x^2+\tilde y^2}) so p \simeq P_0 + Log (2r/(c\varepsilon )). Of course, this corresponds to the expected solution of a sink far enough.

    Now we write the asymptotic matching \displaystyle lim_{\tilde y \rightarrow \infty} \tilde p =lim_{y \rightarrow 0} p

    this gives : P_0 + Log(2) -Log( \varepsilon)= Log(\pi) which provides the value of pressure P_0 that we have to put at the small slit

    Velocity

    Velocity is constant and is \pi at the top of the silo. Remember \tilde x=c cosh(\tilde p) cos (\tilde \psi) and \tilde y=c sinh(\tilde p) sin (\tilde \psi) with c=1.

    At the slit, where \psi goes from 0 to \pi and \tilde y=0 and pressure is 0, we have d \tilde x= sin \psi d \psi and velocity is computed with \frac{\partial \tilde y}{\partial \tilde p} so that \frac{\partial \tilde p}{\partial \tilde y} = \frac{1}{sin \psi}. So that the integral is \displaystyle D = \int_{\tilde x =-1}^{\tilde x =1} \frac{\partial \tilde p}{\partial \tilde y} dx = \pi.

    Flow rate in this device is \pi.

    Code

    #include "run.h"
#include "poisson.h"
#define MAXLEVEL 8

scalar p[], source[];
double eps=.1;  // OK Up to 0.025 for level 8
face vector beta[];
mgstats mgp;

    boundary conditions

    far away Log ( sin ( \pi z) ) \simeq Log (exp( i \pi z)/2) = -\pi (i x - y) -Log[2] so that the pressure is indeed linear (\pi y - Log 2 ). We put at the top the exact value 2.44658 indeed note that \pi - Log[2.]= 2.448

    close to the sink we have the overlapping between the two descriptions which is the source so that as \varepsilon=.1

    Log[.1] + Log[\pi] - Log[2]= -1.851 gives the value of pressure that we put at the small slit

    p[right] = neumann(0); 
p[left]  = neumann(0);   
p[top]   = dirichlet(2.4465800) ;  
p[bottom] = fabs(x)<= eps ? dirichlet(log(eps) + log(pi) - log(2)): neumann(0);

    domain is unit

    int main()
{    
    L0=1.;
    Y0=0;
    X0=-L0/2.;
    init_grid (1 << MAXLEVEL);
    run();
}

    coefficient of porosity is constant

    event init (i = 0) {
    foreach_face() {
        beta.x[] = 1; 
    }  
}

    no source

    event defaults (i = 0)
{ 
  foreach()
    p[] = source[] = 0.;
  boundary ({p});
}

    At every timestep, but after all the other events for this timestep have been processed (the ‘last’ keyword), we update the pressure field p by solving the Poisson equation with coefficient \beta.

    event pressure (i++, last)
{

    solve \displaystyle \nabla \cdot (\beta \nabla p )= s with http://basilisk.fr/src/poisson.h

      mgp = poisson (p, source, beta);
}

    error

    event logfile (i++)
{
    stats s = statsf (p);
    fprintf (stderr, "%d %g %d %g %g %g\n",
             i, t, mgp.i, s.sum, s.min, s.max);
}

    Save in a file

    event sauve (i++,last)
{
FILE *  fpc = fopen("pressure.txt", "w");
output_field ({p}, fpc, linear = true);
fclose(fpc);
fprintf(stdout," end\n");
}

    Results

    Run

    To compile and run:

     qcc -O2 -Wall -o darcysilo darcysilo.c -lm
./darcyLambSneddon

    or more clean

     make darcysilo.tst; make darcysilo/plots ; make darcysilo.c.html

    Plots

    Results, figures of iso pressure

    set pm3d map
set palette rgbformulae 22,13,-31;
unset colorbox
set xlabel "x  iso p"
 set size (.5*1.6),1
splot [][:] 'pressure.txt' u 1:2:3   not
reset
    (script)

    (script)

    Figure of iso pressure, we see the transition from Lamb solution to constant pressure gradient solution via the infinite sum of sinks from Paterson’s solution.

    L0=1.
reset
set view map
set size (.5*1.6),1
unset key
unset surface
 set contour base
 set cntrparam levels incremental -1.8,.1,2.5
 splot [][:] 'pressure.txt' u 1:2:3  w l not
    (script)

    (script)

    Along x=0 we have for the outer solution Re[Log[i sinh(\pi y)) for the complex potential, so that it is for pressure log(sinh(\pi y)). The linear solution is \pi y-log(2) (in black), approximation for the outer solution. We plot (in red) the numerical solution along x=0 compared with analytical outer solution log(sinh(\pi y)) (in green) for a sink, they diverge near the slit as expected. Lamb inner solution (in blue) rescaled by \varepsilon gives the solution near the slit. There is an overlap region for y=O(\varepsilon) where the three are superposed (as expected by the Matched Asymptotic Expansion theory).

    set key bottom
set xlabel "y"
set ylabel "p(0,y)"

# plot [:][:] 'pressure.txt' u 2:(abs($1)<.01?($3):NaN)  t'num.',log(sinh(pi*x)) 
 plot [:1][:] 'pressure.txt' u 2:(abs($1)<.01?($3):NaN)  t'num.',log(sinh(pi*x)),\
 '../darcyLambSneddon/pressure.txt' u ($2*.1):(abs($1)<.01?($3-1.851):NaN)  t'num. Lamb' w l,pi*(x)-log(2) t'lin. approx' w l lc -1
    (script)

    (script)

    Composite expansion: \displaystyle p_{composite} = p_{inner}+p_{outer} - p_{common\;limit}

    along the center line x=0, we have p_{inner}=\tilde p= log(\pi)-log(2)+log(\varepsilon)+asinh(y/\varepsilon), and p_{outer}=log(sh(\pi y)), and the common behavior that we remove is the sink log(y)+log(\pi) so that we plot p_{composite} compared to the numerics:

    set key bottom
set xlabel "y"
set ylabel "p(0,y)"
eps=.1
plot [:1][:] 'pressure.txt' u 2:(abs($1)<.01?($3):NaN) t'num.'\
,log(pi)-log(2)+log(eps)+asinh(x/eps)+log(sinh(pi*x))-log(x)-log(pi) t'composite'
    (script)

    (script)

    Bibliography