Newton method for the stationary Navier-Stokes equations doesn't converge

4 months ago by
Dear all,

My question is, what would be an optimal approach while trying to solve stationary Navier-Stokes equations in the incompressible case.

So far, I tried to implement my problem using the NonlinearVariationalSolver. In the case of Stokes equations it converges within one iteration, just as it should, however in the Navier-Stokes case it clearly doesn't converge. This probably indicates that I should do it differently. I tried to look for examples within tutorials. Unfortunately everything that i could find was about the non-stationary case, mostly using the Chorin's projection method. I'm not sure whether i could also apply it in my case.

Another example are the Stokes equations here:

Again, is it possible to apply something similar to the Navier-Stokes problem also?

Here is a piece of code that i wrote. Firstly solve the Stokes problem which works fine, then pass the solution to the NonlinearVariationalSolver again as the initial guess for the Navier-Stokes problem.

        def navier_stokes(NS, U_):

                # Define variational problem
                U = TrialFunction(V)
                (u, p) = split(U)
                V_ = TestFunction(V)
                (v, q) = split(V_) 
                (u_, p_) = split(U_)
                F = A(u, v) + B(p, v) + bool(NS)*C(u, u, v) + B(q, u) - L(v)
                F = action(F, U_)
                J = derivative(F, U_, U)
                problem = NonlinearVariationalProblem(F, U_, bcs, J)
                # Create Newton solver
                solver = NonlinearVariationalSolver(problem)
                return U_

        # Call navier-stokes and get sub-functions
        U_ = Function(V)
        U_s = navier_stokes(0, U_)
        (u_s, p_s) = U_s.split(U_s)
        U_ns = navier_stokes(1, U_s)
        (u_ns, p_ns) = U_ns.split(U_ns)​
The variational forms are very classical

\[ \begin{cases}
A(\textbf{u}, \textbf{v}): = \nu \int_{\Omega}\nabla \textbf{u}\cdot \nabla\textbf{v}\hspace{0.1 cm}d\Omega &\textbf{u},\textbf{v}\in V_1\\
B(\textbf{u}, q): =-\frac{1}{\rho} \int_{\Omega}div(\textbf{u})q \hspace{0.1 cm}d\Omega & \textbf{u} \in V_1, q \in V_2 \\
C(\textbf{w};\textbf{z},\textbf{v}): = \int_{\Omega} \sum_{i,j=1}^{d} w_j \frac{\partial z_i}{\partial x_j}v_i \hspace{0.1 cm}d\Omega & \textbf{w},\textbf{z},\textbf{v} \in V_1 \\
L(\textbf{v}): = \frac{1}{\rho}\int_{\Omega}\textbf{f}\cdot \textbf{v}\hspace{0.1 cm}d\Omega &\textbf{v} \in V_1
\end{cases} \]

Here is the picture of the Stokes problem

When it comes to the boundary conditions - there are two parabolic inflows on the left, do-nothing condition on the right and noslip condition everywhere else.

I use classical Taylor-Hood elements.

Community: FEniCS Project
Check out Jan Blechta's tutorial:
written 4 months ago by Adam Janecka  
Thank you so much for your comment! However i am already familiar with this tutorial and they simply use the command solve(F == 0, w, bcs) for the variational form and unfortunately this doesn't work; I tried already. Maybe because my domain is a bit more complicated.
written 4 months ago by Martyna Soszyńska  
Here's a working code for the steady incompressible Navier--Stokes Poiseuille flow in 2D channel using the NonlinearVariationalSolver. Hopefully, it can help.
from __future__ import print_function, division
from fenics import *

parameters['form_compiler']['representation'] = 'uflacs'
parameters['form_compiler']['optimize'] = True

mu = Constant(0.5) # Viscosity
C  = 1.0e4 # Pressure gradient

N = 16
mesh = RectangleMesh(Point(0, -1), Point(10, 1), 5*N, N, 'crossed')

inflw_bndr   = 'near(x[0], 0.0)'
outflw_bndr  = 'near(x[0], 10.0)'
no_slip_bndr = 'near(x[1], -1) || near(x[1], 1)'

### Space definition
V  = VectorElement('P', mesh.ufl_cell(), 2)
P  = FiniteElement('P', mesh.ufl_cell(), 1)
TH = MixedElement([V, P])
W  = FunctionSpace(mesh, TH)

### Define Boundary Conditions
exact_vlct = Expression(('(C/2.0/mu)*(1.0-pow(x[1], 2.0))', '0.0'), C=C, mu=mu, degree=2)

no_slip_bc  = DirichletBC(W.sub(0), Constant((0.0, 0.0)), no_slip_bndr)
u_inflw_bc  = DirichletBC(W.sub(0), exact_vlct, inflw_bndr)
p_outflw_bc = DirichletBC(W.sub(1), Constant(0.0), outflw_bndr)

bcs = [no_slip_bc, u_inflw_bc, p_outflw_bc]

### Steady Part of the Momentum Equation
def steady(u):
    T = -p*I + 2*mu*sym(grad(u))
    return ( inner(grad(u)*u, v_) + inner(T, grad(v_)) - inner(f, v_) ) * dx

### Unknown and test functions
(v_, p_) = TestFunctions(W)
w = Function(W)
(v, p) = split(w)

f = Constant((0.0, 0.0))

I = Identity(v.geometric_dimension())

F = steady(v) + p_*div(v)*dx

J = derivative(F, w)
problem = NonlinearVariationalProblem(F, w, bcs, J)
solver  = NonlinearVariationalSolver(problem)

### Create Files for Storing the Solution
vfile = XDMFFile('results/2d-stationary-channel-ns-velocity.xdmf')
pfile = XDMFFile('results/2d-stationary-channel-ns-pressure.xdmf')
### Compute Solution

### Extract Solutions
(v, p) = w.split()
v.rename('v', 'velocity')
p.rename('p', 'pressure')

written 4 months ago by Adam Janecka  

1 Answer

4 months ago by

solving a fluid dynamics problem of an incompressible Newtonian fluid modeled by the Navier--Stokes material equation is a big challenge. There are different strategies like staggering scheme (Chorin's projection), limiting procedure (like stabilized FEM), or a direct numerical solution (DNS). Often you need to fulfill the inf-sup condition (by using Taylor--Hood elements) as well.

If you are interested for something else, have a look here:
The code is written in FEniCS and you can get it here:

In case you have a question about the method or implementation, if it works or doesn't work, write me an e-mail.

Best, Emek
Please login to add an answer/comment or follow this question.

Similar posts:
Search »