2D buoyancy-driven flow (Boussinesq) through inlet patch at the top wall with circular obstacle in the room

Errors
1

from fenics import *
import numpy as np

-------------------------

Parameters

-------------------------

rho = 1.0 # density
mu = 0.01 # viscosity
kappa = 0.01 # thermal diffusivity
beta = 0.01 # thermal expansion
g = Constant((0.0, -9.81)) # gravity
T_in = 290.0
T_hot = 310.0
dt = 0.01
t_end = 0.5

-------------------------

Mesh

-------------------------

nx, ny = 40, 40
mesh = UnitSquareMesh(nx, ny)

-------------------------

Function spaces

-------------------------

V = VectorElement(“P”, mesh.ufl_cell(), 2) # velocity
Q = FiniteElement(“P”, mesh.ufl_cell(), 1) # pressure
S = FiniteElement(“P”, mesh.ufl_cell(), 1) # temperature
W = FunctionSpace(mesh, MixedElement([V, Q, S]))

w = Function(W) # current solution
w_n = Function(W) # previous solution
(u, p, T) = split(w)
(u_n, p_n, T_n) = split(w_n)
(v, q, s) = TestFunctions(W)

-------------------------

Boundary conditions

-------------------------

inlet = “near(x[0], 0.0)”
outlet = “near(x[0], 1.0)”
walls = “near(x[1], 0.0) || near(x[1], 1.0)”

bcu_inlet = DirichletBC(W.sub(0), Constant((1.0, 0.0)), inlet)
bcu_walls = DirichletBC(W.sub(0), Constant((0.0, 0.0)), walls)
bcp_outlet = DirichletBC(W.sub(1), Constant(0.0), outlet)
bcT_inlet = DirichletBC(W.sub(2), Constant(T_in), inlet)
bcT_hotwall = DirichletBC(W.sub(2), Constant(T_hot), “near(x[1], 0.0)”)
bcs = [bcu_inlet, bcu_walls, bcp_outlet, bcT_inlet, bcT_hotwall]

-------------------------

Weak forms

-------------------------

Time discretization

u_mid = 0.5*(u + u_n)
T_mid = 0.5*(T + T_n)

Momentum

F_m = rhodot((u - u_n)/dt, v)dx
+ rhodot(dot(u_mid, nabla_grad(u_mid)), v)dx
+ 2muinner(sym(grad(u_mid)), sym(grad(v)))dx
- div(v)pdx
- qdiv(u)dx
- rhobeta*(T_mid - T_in)*dot(g, v)*dx

Energy

F_T = (T - T_n)/dtsdx
+ dot(u_mid, grad(T_mid))sdx
+ kappa*dot(grad(T_mid), grad(s))*dx

F = F_m + F_T

-------------------------

Stabilization (SUPG/PSPG)

-------------------------

h = CellDiameter(mesh)
u_mag = sqrt(dot(u_mid, u_mid) + 1e-10)
tau_m = h/(2.0u_mag) # momentum stabilization
tau_c = hh/(2.0*kappa) # thermal stabilization

r_m = rho*((u - u_n)/dt + dot(u_mid, nabla_grad(u_mid)))
- div(2musym(grad(u_mid))) + grad(p)
- rhobeta(T_mid - T_in)g
r_T = (T - T_n)/dt + dot(u_mid, grad(T_mid)) - kappadiv(grad(T_mid))

F += tau_m*dot(r_m, dot(u_mid, grad(v)))*dx \

  • tau_cr_Tdot(u_mid, grad(s))*dx

-------------------------

Solve nonlinear system

-------------------------

J = derivative(F, w)

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

prm = solver.parameters
prm[“nonlinear_solver”] = “newton”
prm[“newton_solver”][“absolute_tolerance”] = 1e-8
prm[“newton_solver”][“relative_tolerance”] = 1e-7
prm[“newton_solver”][“maximum_iterations”] = 25
prm[“newton_solver”][“linear_solver”] = “mumps”

-------------------------

Time stepping

-------------------------

t = 0.0
vtkfile_u = File(“velocity.pvd”)
vtkfile_T = File(“temperature.pvd”)

while t < t_end:
t += dt
print(f"Time = {t:.3f}")
solver.solve()
(u_sol, p_sol, T_sol) = w.split()
vtkfile_u << (u_sol, t)
vtkfile_T << (T_sol, t)
w_n.assign(w)

My solver keeps getting diverging

2

Please format you code appropriately with

```python
# Add code here
```

and also add the output of your program as

```bash
$ Add output here
```
3

Python

from fenics import *
from mshr import *
import numpy as np
import os

# --------------------------
# Simulation parameters
# --------------------------
rho = 1.0
mu = 0.01
U_in = 0.01
dt = 0.001
T = 0.5
num_steps = int(T/dt)
ramp_time = 0.2

# Obstacle
obstacle_center = (0.5, 0.5)
obstacle_radius = 0.1

# Mesh refinement
base_res = 64
refine_radius = 0.15

# Arc inlet parameters
theta_start = 295    # degrees
theta_end = 360  # degrees
inlet_thickness = 0.01

# --------------------------
# Create geometry with mshr
# --------------------------
domain = Rectangle(Point(0.0, 0.0), Point(1.0, 1.0)) - Circle(Point(*obstacle_center), obstacle_radius, 64)
mesh = generate_mesh(domain, base_res)

# --------------------------
# Refine mesh near obstacle
# --------------------------
cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
cell_markers.set_all(False)
for cell in cells(mesh):
    x, y = cell.midpoint().x(), cell.midpoint().y()
    if (x - obstacle_center[0])**2 + (y - obstacle_center[1])**2 < refine_radius**2:
        cell_markers[cell] = True
mesh = refine(mesh, cell_markers)

# --------------------------
# Mixed function space
# --------------------------
V_elem = VectorElement("P", mesh.ufl_cell(), 2)
Q_elem = FiniteElement("P", mesh.ufl_cell(), 1)
W_elem = MixedElement([V_elem, Q_elem])
W = FunctionSpace(mesh, W_elem)
V_space = W.sub(0).collapse()
Q_space = W.sub(1).collapse()

# --------------------------
# ParaView output
# --------------------------
os.makedirs("cavity_flow", exist_ok=True)
vtkfile_u = File("cavity_flow/velocity.pvd")
vtkfile_p = File("cavity_flow/pressure.pvd")

# --------------------------
# Boundary definitions
# --------------------------
tol = 1E-14

class BottomWall(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1],0.0, tol)

class LeftWall(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],0.0, tol)

class RightWall(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],1.0, tol)

bottom_wall = BottomWall()
left_wall = LeftWall()
right_wall = RightWall()

# --------------------------
# Arc-shaped radial inlet
# --------------------------
class RadialArcInlet(UserExpression):
    def __init__(self, center, r_in, theta_start, theta_end, U_in, **kwargs):
        super().__init__(**kwargs)
        self.xc, self.yc = center
        self.r_in = r_in
        self.theta_start = np.deg2rad(theta_start)
        self.theta_end = np.deg2rad(theta_end)
        self.U_in = U_in

    def eval(self, values, x):
        dx = x[0] - self.xc
        dy = x[1] - self.yc
        r = np.sqrt(dx**2 + dy**2)
        if r == 0.0:
            values[0] = values[1] = 0.0
            return
        theta = np.arctan2(dy, dx)
        # Normalize theta to 0-2pi
        if theta < 0.0:
            theta += 2*np.pi
        if self.theta_start <= theta <= self.theta_end and abs(r - self.r_in) <= inlet_thickness:
            values[0] = self.U_in * dx / r
            values[1] = self.U_in * dy / r
        else:
            values[0] = 0.0
            values[1] = 0.0

    def value_shape(self):
        return (2,)

# Create inlet expression
inlet_expr = RadialArcInlet(obstacle_center, obstacle_radius, theta_start, theta_end, U_in, degree=2)

# --------------------------
# Boundary conditions
# --------------------------
bcu_inlet = DirichletBC(W.sub(0), inlet_expr, 'on_boundary')  # Apply to entire boundary; velocity = 0 elsewhere automatically
bcu_bottom = DirichletBC(W.sub(0), Constant((0,0)), bottom_wall)
bcu_left = DirichletBC(W.sub(0).sub(0), Constant(0.0), left_wall)
bcu_right = DirichletBC(W.sub(0).sub(0), Constant(0.0), right_wall)
bcs = [bcu_inlet, bcu_bottom, bcu_left, bcu_right]

# --------------------------
# Functions
# --------------------------
w = Function(W)
w_n = Function(W)
(u, p) = split(w)
(u_n, p_n) = split(w_n)
(v, q) = TestFunctions(W)

# Initialize velocity + zero pressure
u_init = interpolate(inlet_expr, V_space)
p_init = interpolate(Constant(0.0), Q_space)
assign(w_n, [u_init, p_init])

# --------------------------
# Variational form (BDF1)
# --------------------------
F = (rho*dot((u - u_n)/dt, v)*dx
     + rho*dot(dot(u_n, nabla_grad(u)), v)*dx
     + mu*inner(grad(u), grad(v))*dx
     - div(v)*p*dx
     + q*div(u)*dx)

# --------------------------
# Time-stepping
# --------------------------
t = 0.0
for n in range(num_steps):
    t += dt
    ramp_factor = min(1.0, t/ramp_time)
    inlet_expr.U_in = U_in * ramp_factor

    # First step linearization for stability
    if n == 0:
        F_linear = (rho*dot((u - u_n)/dt, v)*dx
                    + mu*inner(grad(u), grad(v))*dx
                    - div(v)*p*dx
                    + q*div(u)*dx)
        problem = NonlinearVariationalProblem(F_linear, w, bcs, J=derivative(F_linear, w))
    else:
        problem = NonlinearVariationalProblem(F, w, bcs, J=derivative(F, w))

    solver = NonlinearVariationalSolver(problem)
    solver.parameters['newton_solver']['absolute_tolerance'] = 1E-8
    solver.parameters['newton_solver']['relative_tolerance'] = 1E-6
    solver.parameters['newton_solver']['maximum_iterations'] = 50
    solver.parameters['newton_solver']['report'] = True
    solver.parameters['newton_solver']['linear_solver'] = 'mumps'

    solver.solve()

    (u_sol, p_sol) = w.split()
    vtkfile_u << (u_sol, t)
    vtkfile_p << (p_sol, t)
    w_n.assign(w)

print("Simulation completed.")

Newton iteration 49: r (abs) = 5.000e-03 (tol = 1.000e-08) r (rel) = 3.768e+00 (tol = 1.000e-06)
Newton iteration 50: r (abs) = 5.129e-03 (tol = 1.000e-08) r (rel) = 3.864e+00 (tol = 1.000e-06)
Traceback (most recent call last):
File “CavityFlow5.py”, line 186, in
solver.solve()
RuntimeError:

*** -------------------------------------------------------------------------

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0

*** DOLFIN version: 2019.1.0
*** Git changeset: 828b3965c3ebcd128da494c596c78b17c0a303da
*** -------------------------------------------------------------------------

4

A problem here is that you use mshr, which is not available on most systems any more.
Have you considered switching to DOLFINx? As legacy FEniCS has not been maintained for over 5 years now.

Secondly, as your problem runs for a while, have you looked at the solution and seen if you can spot an instability building up in the proximity of the circular object?