This is a viscous flow described with the Navier-Stokes equation. The top is a free surface that can deform accordingly to the pressure jump through the interface due to surface tension, and the normal and tangential stress in the flow. At the inflow we impose a profile that looks like a boundary layer.

    clear all; clf; format compact
    for tre=1:length(paramvec);
        % parameters
        Re=500; % Reynolds number
        Nx=100;  % number of grid nodes in z
        Ny=50;  %number of grid nodes in r
        Lx=10;  % domain length
        Ly=1;   % domain height
        sig=5;  % surface tension
        V=Lx*Ly*param; % desired volume
        umax=1; % max inflow velocity
        % differentiation
        % useful[l.ctl;; l.ctr]; % include the corners to the top
        l.u=(1:NN)'; l.v=l.u+NN; l.p=l.v+NN; l.e=l.p(end)+(1:Nx)'; l.p0=l.e(end)+1; % location vectors
        l.ue=[l.u;l.e];[l.v;l.e];[l.p;l.e]; % compound location vectors
        p0loc=2*Ny; % where to impose pressure=p0; % a shortcut because its often used
        n=3*NN+Nx+1; % numbers of degrees of freedom
        ZZ=spalloc(n,n,0);    zx=spalloc(Nx,Nx,0);
        II=speye(n);        ix=speye(Nx);
        Iu=II(l.u,:); Iv=II(l.v,:); Ip=II(l.p,:); Ie=II(l.e,:); Ip0=II(l.p0,:);
        T=kron(speye(Nx,Nx),ones(Ny,1)); % transformation operator from 1D x to 2D grid
        % initial guess
        p=-mu*X*umax; p=p-p(p0loc);
        % Newton iterations
        disp('Newton loop'); tlu=0;
        if tre==1; q=q0; end
        while ~quit
            % the present solution and its computational-space derivatives
            e=q(l.e); ex=d.x*e; exx=d.xx*e; E=T*e; Ex=T*ex; Exx=T*exx; a=1+ex.^2;
            u=q(l.u); v=q(l.v); p=q(l.p); p0=q(l.p0);
            ux=Dc.x*u; uy=Dc.y*u; uxy=Dc.y*ux; uyy=Dc.yy*u;
            vx=Dc.x*v; vy=Dc.y*v; vxy=Dc.y*vx; vyy=Dc.yy*v;
            px=Dc.x*p; py=Dc.y*p;

    Physical space differentiation matrices

    This is just like we did in stretching_formula.m.

            % physical-space differentation matrices
            D.lap=Dc.xx ...
                +spd(-2*yy.*E.^-1.*Ex)*(Dc.x*Dc.y) ...
                +spd(yy.^2.*E.^-2.*Ex.^2+E.^-2)*Dc.yy ...

    Jacobians of the differentiation matrices

    Since the domain deformation is one of the unknown of th system, the differentiation is no longer a linear operation (because it depends nonlinearly on the domain stretching \eta). We do just like we did in the test case jacobian_formula.m.

            % jacobians of the physical-space differentiation operators
            D.xue=spd(yy.*uy)*DD;  D.yue=spd(uy)*EE;  D.lapue=spd(uxy)*AA+spd(uyy)*BB+spd(uy)*CC;
            D.xve=spd(yy.*vy)*DD;  D.yve=spd(vy)*EE;  D.lapve=spd(vxy)*AA+spd(vyy)*BB+spd(vy)*CC;
            D.xpe=spd(yy.*py)*DD;  D.ype=spd(py)*EE;
            Iux=[D.x,D.xue]*II(l.ue,:);  Iuy=[D.y,D.yue]*II(l.ue,:);  Iulap=[D.lap,D.lapue]*II(l.ue,:);
            Ivx=[D.x,D.xve]*II(,:);  Ivy=[D.y,D.yve]*II(,:);  Ivlap=[D.lap,D.lapve]*II(,:);
            Ipx=[D.x,D.xpe]*II(,:);  Ipy=[D.y,D.ype]*II(,:);
            % physical-space derivatives
            ux=D.x*u; uy=D.y*u;
            vx=D.x*v; vy=D.y*v;
            px=D.x*p; py=D.y*p;

    The nonlinear function

    The first thing we do is to compute the expression of the normal and tangential stresses at the free surface. The stress tensor in a viscous Newtonian fluid is \displaystyle \sigma= \begin{pmatrix} 2\mu u_x-p & \mu(u_y+v_x) \\ \mu(u_y+v_x) & 2\mu v_y-p \end{pmatrix} We are interested at the components of the stress along our free-surface. We first need to find the expresion of its normal and tangent vectors. We can say that the normal is along the gradient of a function for which the free-surface is an iso-line. For instance take \displaystyle F(x,y)=y-\eta(x) then its gradient, normalized to unit one is the normal to the interface \displaystyle \overrightarrow{n}= \begin{pmatrix} -\frac{\eta_x}{\sqrt{1+\eta_x^2}}\\ \frac{1}{\sqrt{1+\eta_x^2}}\\ \end{pmatrix} and the tangent vector is \displaystyle \overrightarrow{t}= \begin{pmatrix} \frac{1}{\sqrt{1+\eta_x^2}}\\ \frac{\eta_x}{\sqrt{1+\eta_x^2}}\\ \end{pmatrix} Then, the normal stress is \displaystyle \sigma_n=\overrightarrow{n}\sigma \overrightarrow{n} =\frac{\eta_x^2(2\mu u_x-p)-2\mu\eta_x(u_y+v_x)+2\mu v_y-p}{1+\eta_x^2} and the tangential stress is \displaystyle \sigma_t=\overrightarrow{t}\sigma \overrightarrow{n} =\mu\frac{(1-\eta_x^2)(u_y+v_x)+2\mu\eta_x(v_y-u_x)}{1+\eta_x^2} The conditions at the interface are that the normal stress has a jump equal to the pressure jump through a curved interface with surface tension, and the tangential stress is zero.

    We have then Navier-Stokes in the flow domain. And the last equation imposes the value of the fluid volume.

            % nonlinear function
            stressn=-sig*(exx.*a.^-1.5)+(ex.^2.*(2*mu*ux(t)-p(t))-2*mu*ex.*(uy(t)+vx(t))+2*mu*vy(t)-p(t))./a; % normal stress
            stresst=mu*((1-ex.^2).*(uy(t)+vx(t))+2*ex.*(vy(t)-ux(t)))./a; % tangential stress
            f=[-(u.*ux+v.*uy)+D.lap*u/Re-px; ...
                -(u.*vx+v.*vy)+D.lap*v/Re-py; ...
                ux+vy; ...
                stressn; ... % jump in normal stress
                d.wx*Ly*e-V]; % volume

    The Jacobian

    We build here the expression of the Jacobian of f.

            % analytical Jacobian
            Astressn=-sig*(spd(a.^-1.5)*d.xx-spd(3*ex.*exx.*a.^-2.5)*d.x)*Ie ...
                +diag(2*mu*ex.^2./a)*Iux(t,:)-diag(2*mu*ex./a)*Iuy(t,:) ...
                -diag(2*mu*ex./a)*Ivx(t,:)+diag(2*mu./a)*Ivy(t,:) ...
            Astresst=mu*(diag((1-ex.^2)./a)*Iuy(t,:)-diag(2*ex./a)*Iux(t,:) ...
            A=[-(spd(u)*Iux+spd(ux)*Iu+spd(v)*Iuy+spd(uy)*Iv)+Iulap/Re-Ipx; ...
                -(spd(u)*Ivx+spd(vx)*Iu+spd(v)*Ivy+spd(vy)*Iv)+Ivlap/Re-Ipy; ...
                Iux+Ivy; ...
                Astressn; ...

    Boundary conditions

    They are:

    • Dirichlet for u and v at inflow (imposed inflow profile set by the initial guess) and bottom (no-slip)
    • Neumann for u and v at outflow (x derivative is zero)
    • The pressure must equal the reference pressure p_0 close to the bottom left corner
    • The Neumann condition on pressure to avoid spurious pressure modes
    • The height of the interface is imposed at x=0 and x=L_x.
    • The free surface is a streamline, or said differently: no penetration through the interface: \displaystyle v-u\eta_x=0
    • The tangential stress is zero
            % Boundary conditions
            dir=[l.u([l.left;]); l.v([l.left;])]; % where to put Dirichlet on u and v
            neuploc=[l.ctl;l.ctr;l.ctr-Ny];  % where to impose the neumann condition on the pressure
            loc=[dir; ...
                l.u(l.right); ...
                l.v(l.right); ...
                l.p(p0loc); ...
                l.p(neuploc); ...
                l.e([1,Nx]); ...
            f(loc)=[q(dir)-q0(dir); ... % all Dirichlet
                ux(l.right); ... % Homogeneous Neuman on u at outflow
                vx(l.right); ...
                p(p0loc)-p0; % pressure equal to p0 at p0loc
                D.lap(neuploc,:)*v/Re-px(neuploc); ... % neumann constraint on pressure
                e(1)-Ly; ... % free-surface pinned at both ends
                e(Nx)-Ly; ...
                v(t)-u(t).*ex; ... % free-surface is a streamline
            C=[II(dir,:); ...     % all Dirichlet
                Iux(l.right,:); ...
                Ivx(l.right,:); ...
                Ip(p0loc,:)-Ip0;  % pressure equal to p0 at p0loc
                Ivlap(neuploc,:)/Re-Ipx(neuploc,:); ... % Neumann constraint on pressure
                ix(1,:)*Ie; ... % free-surface pinned at both ends
                ix(Nx,:)*Ie; ...
                Iv(t,:)-diag(ex)*Iu(t,:)-spd(u(t))*d.x*Ie; ... % free-surface is a streamline
            % plotting
            Yp=Y.*reshape(E,Ny,Nx); % physical mesh
            ww=reshape(u,Ny,Nx); vv=reshape(v,Ny,Nx);
            vo=vx-uy; % vorticity
            surf(X,Yp,reshape(vo,Ny,Nx)-10,'facealpha',1); view(2); shading interp; hold on
            sely=1:2:Ny; selx=round(linspace(1,Nx,10));
            xlabel('x'); ylabel('y'); title('Vorticity and velocity field'); grid off;hold off
            % convergence test
            disp([num2str(count) '  ' sprintf('%9.3g',res) '  tjac:' sprintf('%7.3g',tjac) '  tlu:' sprintf('%7.3g',tlu)]);
            if count>50||res>1e5; disp('no convergence');break; end
            if res<1e-5; quit=1; disp('converged'); continue; end
            % Newton step
            tic; q=q-A\f; tlu=toc;

    Here is the final figure:

    The figure

    The figure