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))


    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 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.


    #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()
        init_grid (1 << MAXLEVEL);

    coefficient of porosity is constant

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

    no source

    event defaults (i = 0)
        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

      mgp = poisson (p, source, beta);


    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);
    fprintf(stdout," end\n");



    To compile and run:

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

    or more clean

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


    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


    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.

    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


    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


    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)"
    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'