# Uzawa scheme for stokes bingham flow

**URL:** <https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233>\
**Category:** General Discussion\
**Created:** [May 9, 2024, 11:10am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233 "2024-05-09T11:10:13Z")\
**Posts on this page:** 1\
**Showing post:** 16

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [May 13, 2024, 1:15pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/16 "2024-05-13T13:15:19Z")

</div>

In that paper the flow is taken incompressible, whereas in your code you use compressible. This is why the solution is different (but looks correct).  
You can modify easily your code to do the incompressible case, by using the extra variables p,q as  
`solve bing([ux,uy,p],[vx,vy,q])`  
and adding  
`-int2d(Th)(p*Divergence(vx, vy)+q*Divergence(ux, uy))`  
But then it happens that the scheme converges extremely slowly (you need thousands of iterations). In order to improve it you should use the Augmented Lagrangian method as in that paper.

* * *

Another point is that your computation of the stream function could be misleading. You look for \psi satisfying  
-\Delta\psi=\partial\_y u\_x-\partial\_x u\_y.  
An important property that one would like that if \partial\_x u\_x+\partial\_y u\_y=0 then \partial\_x\psi=u\_y and -\partial\_y\psi=u\_x.  
You see then that the Dirichlet condition on \psi is not appropriate. It is better to take  
\frac{\partial\psi}{\partial n}=u\_y n\_x-u\_x n\_y on \partial\Omega.  
We are thus led to the problem  
-\Delta \psi=f in \Omega,\qquad\frac{\partial\psi}{\partial n}=g on \partial\Omega.  
To solve this system you need the compatibility condition (to see that it is necessary, integrate the equation over \Omega)  
\int\_\Omega f=-\int\_{\partial\Omega}g  
Assuming this, the system is not invertible since \psi is defined up to a constant. You can determine it by imposing \int\_{\partial\Omega}\psi=0.  
In order to get the solution you can then solve  
-\Delta \psi=f in \Omega,\qquad\frac{\partial\psi}{\partial n}+\varepsilon\psi=g on \partial\Omega,  
for some small \varepsilon. This system is well-posed. This leads to the code

```auto
fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
    = int2d(Th)(
          dx(dyg)*dx(dyg2)
        + dy(dyg)*dy(dyg2)
    )
    + int2d(Th)(
        - dyg2*(dy(ux) - dx(uy))
    )
    + int1d(Th)(
        - dyg2*(uy*N.x-ux*N.y)
    )
    + int1d(Th)(1e-8*dyg*dyg2)
    ;

```

---

_[View the full topic](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233)._