Is the variational formulation correct?

variational formulation
1

The following code is used to solve the incompressible flow. The computaional domain is as follows

图片1
图片11352×392 23.4 KB

My question: the variational formulation of the NS equation is correct?
Thank you very much!

from dolfin import *
import mshr
import matplotlib.pyplot as plt
from mshr import *
from fenics import *

mu = 1.05E-3
rho = 1000
#Generate mesh with mshr

domain = (Rectangle(dolfin.Point(0., 0.), dolfin.Point(15.0E-3, 5.0E-3)) +
Rectangle(dolfin.Point(15.0E-3, 0.), dolfin.Point(135.0E-3, 30.0E-3)) +
Rectangle(dolfin.Point(135.0E-3, 0.), dolfin.Point(150.0E-3, 5.0E-3)))
mesh = generate_mesh(domain, 100)
V = VectorElement('P', mesh.ufl_cell(), 2)
Q = FiniteElement('P', mesh.ufl_cell(), 1)
TH = MixedElement([V, Q])
W = FunctionSpace(mesh, TH)
#Defining the boundary
#plot(mesh)
inflow = 'near(x[0], 0.0)'
outflow = 'near(x[0], 150.0E-3)'
walls = ('near(x[1], 5.0E-3) || near(x[1], 30.0E-3) || near(x[0],15.0E-3) ||' + 'near(x[0], 135.0E-3)')
symm = 'near(x[1], 0.0)'

#Defining boundary conditions
#plot(mesh)
bcu_inflow = DirichletBC(W.sub(0), Constant((0.01, 0.0)), inflow)
bcp_outflow = DirichletBC(W.sub(1), Constant(0), outflow)
bcu_noslip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), walls)
bcu_symm = DirichletBC(W.sub(0).sub(1), Constant(0.0), symm)
bcs = [bcu_noslip, bcu_inflow, bcu_symm, bcp_outflow]
#Define variational forms

v, q = TestFunctions(W)
w = Function(W)
dw= TrialFunction(W)
u, p = split(w)
#Equation

F1 = mu*inner(grad(u),grad(v))*dx + rho*inner(dot(grad(u),u),v)*dx \
     - inner(p, div(v))*dx + inner(q, div(u))*dx 
J = derivative(F1, w, dw)




newton_solver_parameters={"nonlinear_solver": "newton",
                        "newton_solver"     :{"linear_solver"  : "mumps",
                                            "maximum_iterations": 50,
                                            "report": True,
                                            "error_on_nonconvergence": True
                                            }
                        }




problem = NonlinearVariationalProblem(F1, w, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.parameters.update(newton_solver_parameters)

(itera, converged) = solver.solve()


(u, p) = w.split()
plot(u)
2

The form F1 looks like a weak form that is consistent with the Navier–Stokes equations, but you want to be careful about how you integrate by parts in the viscosity term. It is indeed true that, for a solenoidal velocity field \mathbf{u} and spatially-constant viscosity \mu,

-\nabla\cdot(2\mu\operatorname{sym}\nabla\mathbf{u}) = -\mu\Delta\mathbf{u}\text{ ,}

and

\int_\Omega \mu\nabla\mathbf{u} : \nabla\mathbf{v}\,d\Omega

is the interior term most commonly used for the weak Laplacian (for test function \mathbf{v}). However, you want to also consider how the interior terms of the weak form affects Neumann boundary conditions (BCs) in cases where not all components of the velocity are constrained by Dirichlet BCs on the entire boundary (e.g., at outflows, symmetry planes, etc.). Dropping boundary terms from the above integration by parts implies the homogeneous Neumann BC

(\mu\nabla\mathbf{u} - p\mathbf{I})\cdot\mathbf{n} = \mathbf{0}\text{ ,}

where \mathbf{n} is the unit outward normal to \partial\Omega and I’ve assumed the standard integration-by-parts of the pressure term. The physical interpretation of this is somewhat unclear, though. Thus, I typically use the following weak form of the viscosity term:

\int_\Omega (2\mu\operatorname{sym}\nabla\mathbf{u}):\nabla\mathbf{v}\,d\Omega = \int_\Omega 2\mu\operatorname{sym}\nabla\mathbf{u}:\operatorname{sym}\nabla\mathbf{v}\,d\Omega\text{ ,}

where the second equivalent expression emphasizes the form’s symmetry and follows from the applying the identity \operatorname{sym}\mathbf{A}:\operatorname{skew}\mathbf{B} = 0 for rank-2 tensors \mathbf{A} and \mathbf{B} (while the first expression follows more directly from integration by parts). This alternate integration by parts implies the Neumann BC

(2\mu\operatorname{sym}\nabla\mathbf{u} - p\mathbf{I})\cdot\mathbf{n} = \pmb{\sigma}\cdot\mathbf{n} = \mathbf{0}\text{ ,}

where \pmb{\sigma} is the Cauchy stress, i.e., a traction-free BC, which is more physically meaningful.

It’s also generally discouraged to strongly enforce Dirichlet BCs on the pressure, since this breaks the global mass conservation of the weak formulation. (The pressure test space no longer includes a global constant.) Many people simply use traction-free boundary conditions on outflows (which can be implemented by “doing nothing” on the boundary in conjunction with the interior weak form of the viscous term given above), but this can become unstable in cases where there is some possibility of flow reversal on the outflow boundary. There are a variety of techniques for formulating the problem to be robust against this “backflow divergence”, as outlined here, where I typically prefer the method from Section 2.2.1. (Note that this is a modification of the BCs in the continuous problem statement to make it well-posed, not a feature of the discretization.)

Lastly, the Bubnov–Galerkin method implemented in the MWE can become unstable at higher Reynolds numbers, necessitating some form of stabilization for the advection term (e.g., the streamline-upwind Petrov–Galerkin (SUPG) method). If you’re primarily interested in \mu\sim 10^{-3}, \rho\sim 10^3, \vert\mathbf{u}\vert\sim 10^{-2}, and L\sim 10^{-2} (in some consistent system of units), as in the example, this would imply

\text{Re} = \frac{\rho\vert\mathbf{u}\vert L}{\mu} \sim 10^2\text{ ,}

where you can often get away without any stabilization, assuming a moderate mesh resolution. However, it’s worth worrying about if \text{Re} gets into the thousands or higher. (Technically, you want to keep your “element Reynolds number”, where L is roughly the element diameter, low, ideally \lesssim 10^0, but this resolution requirement becomes prohibitive at higher global Reynolds numbers, so stabilization becomes necessary.)

2 Likes
3

Thanks for your suggestion. Thank you! However, is the test function vanished on the symmetric boundary?

4

Yes, your boundary condition is and corresponding weak are correct, regarding the imposition of the symmetry boundary condition. See also my answer to How to impose symmetry conditions on domain centreline, Stokes equations? - #4 by plugged

1 Like
5

However, is the test function vanished on the symmetric boundary?

Only the component of the vector-valued test function \mathbf{v} orthogonal to the symmetry plane goes to zero when strongly enforcing a no-penetration/slip condition (i.e., \mathbf{v}\cdot\mathbf{n} = 0, but \mathbf{v}\neq\mathbf{0}).

1 Like