Using the code from http://math.mit.edu/~gs/cse/codes/mit18086_navierstokes.pdf, I simulated cavity flow with an external force with the boundary conditions on the all 4 sides are $u=v=0$ and $\frac{\partial p}{\partial n}=0$. I just added external force to u direction in the step that calculates nonlinear term $$\text{(1)} \quad \frac{U^*-U^n}{\Delta t} = -((U^n)^2)_x - (U^nV^n)_y $$ to $$\text{(1)} \quad \frac{U^*-U^n}{\Delta t} = -((U^n)^2)_x - (U^nV^n)_y +f $$ . So the code was changed only 2 lines, from

uN = x*0+1
U = U-dt*(UVy(2:end-1,:)+U2x);


uN = x*0
U = U-dt*(UVy(2:end-1,:)+U2x+10);

, and I got this result.enter image description here

I am confused with this result. This external force is like gravity, so I thought, when applying only gravity the fluid inside a cavity does not circulate. It is like water in a glass. Why does it circulate? Is this physically correct or because of the numerical simulation method?

  • $\begingroup$ It might be your pressure boundary conditions. In such a static setup, as you describe it, the gravity -- which can be written as the gradient of a potential field -- is compensated by the pressure gradient. From your setup and from your plot, I would say that $\partial p /\partial n = 0$ is not correct at the uppper and lower boundaries. Can you try to just set $p=0$ at the right (or left) boundary? $\endgroup$
    – Jan
    Commented Nov 9, 2016 at 13:25
  • 1
    $\begingroup$ I can't seem to understand what your formulas represent. Can you state what you choose as a force $f$ when you make these modifications? $\endgroup$ Commented Nov 9, 2016 at 15:41
  • $\begingroup$ I agree with @Jan, your pressure BC does not make sense in the situation you describe. $\endgroup$
    – user20857
    Commented Nov 9, 2016 at 17:48
  • $\begingroup$ @Jan, I changed the boundary condition to set p=0, following the comment "Discretization schemes should not change when you differ the BCs. However the matrices used for multiplication change depending on the boundary conditions, the value at matrix location (1,1) and (n,n)." and the code posted on the link cfd-online.com/Forums/main/90303-matlab-code-pipe-flow.html. But I still get the same result. $\endgroup$
    – user26767
    Commented Nov 9, 2016 at 20:39
  • 1
    $\begingroup$ So what is $f$ then in your example? $\endgroup$ Commented Nov 9, 2016 at 22:25

1 Answer 1


@Wolfgang Bangerth has a very important point about what $f$ is. The constant external force will not affect your solution unless your BCs allow for motion. In your case, they do not. An easy conceptual proof of this is to look at the vorticity equation (by taking the curl of your equations). The source term disappears, which means no motion will ensue. Therefore, the motion you're see is from an non-physical setup.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.