Hello everyone,
I have been trying to implement a turbulence model in FEniCS (I am working in legacy version of FEniCS). So naturally I started with k-\epsilon (k-epsilon) model (see for example here Standard k-epsilon model – CFD-Wiki, the free CFD reference (cfd-online.com)), more precisely the Lam Bremhorst k-\epsilon turbulence model (see here Lam Bremhorst k-ε turbulence model (cham.co.uk)) which differs by using the so called damping functions.
The whole problem (in strong form) is given here:
\frac{\partial u}{\partial t} + (u \cdot \nabla) u = - \nabla p + \nabla \cdot [(\nu + \nu_t) \nabla u]
\nabla \cdot u = 0
\frac{\partial k}{\partial t} + (u \cdot \nabla) k = \nabla \cdot \left[(\nu + \frac{\nu_t}{\sigma_k}) \nabla k \right] + P_k - k R_k
\frac{\partial \epsilon}{\partial t} + (u \cdot \nabla) \epsilon = \nabla \cdot \left[(\nu + \frac{\nu_t}
{\sigma_{\epsilon}}) \nabla \epsilon \right] + P_{\epsilon} - \epsilon R_{\varepsilon}
u = 0, x \in \Gamma_{wall}
p = 20, x \in \Gamma_{in}
p = 0, x \in \Gamma_{out}
k = 0, x \in \Gamma_{wall}
all non-specified b.c. are homogeneous Neumann
Where first two equations just represent the incompressible Navier-Stokes equations with effective viscosity given by \nu_{eff} = \nu + \nu_t. The other two equations are advection-diffusion-reaction equations for two scalar quantities, turbulent kinetic energy (k) and dissipation of turbulent kinetic energy (\epsilon).
The terms that appear in the formulation are as follows:
\nu_t = C_{\nu} f_{\nu} \frac{k^2}{\epsilon}
P_k = \nu_t S^2
R_k = \frac{\epsilon}{k}
P_{\epsilon} = C_1 f_1 R_k P_k
R_{\epsilon} = C_2 f_2 R_k
S = \sqrt{2 \nabla_{s} u:\nabla_{s} u}
Where f_i are the aforementioned damping functions, given by:
f_{\nu} = (1 - \exp(-0.0165 Re_{k}))^2 (1 + \frac{20.5}{Re_t}) *
f_1 = 1 + \left(\frac{0.05}{f_{\nu}} \right)^3
f_2 = 1 - \exp(-Re_t^2)
Re_t = \frac{k^2}{\nu \epsilon}, Re_k = \frac{y\sqrt{k}}{\nu}, y = \text{distance to the wall}
And C_i are the model constants with values
C_{\nu} = 0.09, C_1 = 1.44, C_2 = 1.92
*Note that in my code, function f_{\nu} is implemented such that it is limited from bellow as well as from above.
As you can see most of the terms that appear in my problem are highly nonlinear functions of the problem unknowns (i.e. \nu_t is a function of both k and \epsilon, P_k is also a function of k, \epsilon but also of u and the distance from the wall y and so on for the other terms). That is why I need to linearize my problem such that I use the last known values of u, k, \epsilon (i.e. values from previous iteration) when I calculate those terms (this way I can treat all the equations separately).
This (finally 
) leads me to my question, how should the terms that appear in the variational formulation be implemented? Currently I am doing it like this, which works (for my simple case of turbulent channel flow) but I don’t feel like this is the best way of implementing it.
'''
snippet from my code that is responsible for defining terms that appear
in variational formulation.
- k0, e0 are Functions that are solutions from previous iteration
- u1 is Function that is solution from the current iteration (obtained
from solving NS)
- dfun is an expression that corrsponds to distance to the wall
'''
# Define Min, Max function
def Max(a, b):
return (a+b+abs(a-b))/Constant(2)
def Min(a, b):
return (a+b-abs(a-b))/Constant(2)
# Define helping functions for k-eps model
Re_t = (1. / nu) * (k0**2 / e0)
Re_k = (1. / nu) * (sqrt(k0) * dfun)
f_nu = Max(Min((1 - exp(- 0.0165 * Re_k))**2 * (1 + 20.5 / Re_t), \
Constant(1.0)), Constant(0.01116))
f_1 = 1 + (0.05 / f_nu)**3
f_2 = 1 - exp(- Re_t**2)
S_sq = 2 * inner(sym(nabla_grad(u1)), sym(nabla_grad(u1)))
# Define terms that appear in variational formulation for k-eps model
nut = 0.09 * f_nu * (k0**2 / e0)
prod_k = nut * S_sq
react_k = e0 / k0
prod_e = 1.44 * react_k * f_1 * prod_k
react_e = 1.92 * react_k * f_2
If I understand it correctly, those terms are now UFL terms and one thing that might speed everything would be if they were dolfin functions / dolfin expressions. However when I tried projecting those UFL terms into some FunctionSpace (for example CG1) then my code simulation no longer produced the correct output.
So my second question would be, how to properly project UFL terms onto a FunctionSpace? and how to use such projected functions as, for example, viscosity terms?
** To my knowledge (please correct me if I am wrong), there isn’t much up to date material (with code included) for turbulence modelling in FEniCS so I am posting my whole code for anyone interested in implementing their own turbulence model. This code works and to my knowledge produces the correct results.
# Import dependencies
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
"""
This script simulates the fully developed turbulent channel flow.
Infinitely long channel is simulated using periodic boundary conditions.
Turbulence model that is used is the Lam Bremhorst k-ϵ turbulence model
"""
# Create mesh
height = 2.0; width = 2.0; nx = 4; ny = 64
def meshing(lx,ly, Nx,Ny):
m = UnitSquareMesh(Nx, Ny)
x = m.coordinates()
#Refine on top and bottom sides
x[:, 1] = 0.5 * ( np.arctan(2 * pi * (x[:,1] - 0.5)) / np.arctan(pi) + 1)
#Scale
x[:,0] = x[:,0]*lx
x[:,1] = x[:,1]*ly
return m
mesh = meshing(width, height, nx, ny,)
# Create boundaries subdomains
class Walls(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (near(x[1], 0) or near(x[1], height))
class Inflow(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0)
class Outflow(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], width)
# Class for periodic boundary condition
class Periodic(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0)
# Map right boundary (H) to left boundary (G)
def map(self, x, y):
y[0] = x[0] - width
y[1] = x[1]
# Initialize instances of boundary classes
inflow = Inflow()
outflow = Outflow()
walls = Walls()
periodic = Periodic()
# Choose values for initial and boundary conditions
u_init = 10
k_init = 0.375
e_init = 0.075
# Create function spaces for all variables
V = VectorFunctionSpace(mesh, "CG", 2, constrained_domain = periodic)
Q = FunctionSpace(mesh, "CG", 1)
K = FunctionSpace(mesh, "CG", 1, constrained_domain = periodic)
# Create boundary conditions for u,k,e
bcu = [DirichletBC(V, Constant((0.0, 0.0)), walls)]
bck = [DirichletBC(K, Constant(0.0), walls)]
bce = []
# Create boundary conditions for pressure p
bcp_inflow = DirichletBC(Q, Constant(2.0), inflow)
bcp_outflow = DirichletBC(Q, Constant(0.0), outflow)
bcp = [bcp_inflow, bcp_outflow]
# Define trial functions
u = TrialFunction(V)
p = TrialFunction(Q)
k = TrialFunction(K)
e = TrialFunction(K)
# Define test functions
v = TestFunction(V)
q = TestFunction(Q)
phi = TestFunction(K)
psi = TestFunction(K)
# Define functions for quantities at previous time
u0 = interpolate(Constant((u_init, 0)), V)
k0 = interpolate(Constant(k_init), K)
e0 = interpolate(Constant(e_init), K)
# Define functions where solutions will be stored
u1 = Function(V)
p1 = Function(Q)
k1 = Function(K)
e1 = Function(K)
# Define constantss used in variational forms
dt = Constant(0.05)
nu = Constant(0.00181818)
# Calculate distance to the wall
dfun = Expression('1 - abs(1 - x[1])', degree = 1)
# Define Min, Max function
def Max(a, b):
return (a+b+abs(a-b))/Constant(2)
def Min(a, b):
return (a+b-abs(a-b))/Constant(2)
####################################### MAIN PART #######################################
# Define helping functions for k-eps model
Re_t = (1. / nu) * (k0**2 / e0)
Re_k = (1. / nu) * (sqrt(k0) * dfun)
f_nu = Max(Min((1 - exp(- 0.0165 * Re_k))**2 * (1 + 20.5 / Re_t), \
Constant(1.0)), Constant(0.01116))
f_1 = 1 + (0.05 / f_nu)**3
f_2 = 1 - exp(- Re_t**2)
S_sq = 2 * inner(sym(nabla_grad(u1)), sym(nabla_grad(u1)))
# Define terms that appear in variational formulation for k-eps model
nut = 0.09 * f_nu * (k0**2 / e0)
prod_k = nut * S_sq
react_k = e0 / k0
prod_e = 1.44 * react_k * f_1 * prod_k
react_e = 1.92 * react_k * f_2
#########################################################################################
# Create variational formulation for N-S (Chorin's splitting method)
F1 = (1. / dt) * dot(u - u0, v)*dx \
+ dot(dot(u0, nabla_grad(u)), v)*dx \
+ (nu + nut) * inner(grad(u), grad(v))*dx
a_1 = lhs(F1); l_1 = rhs(F1)
a_2 = dot(grad(p), grad(q))*dx
l_2 = (-1. / dt) * dot(div(u1), q)*dx
a_3 = dot(u, v)*dx
l_3 = dot(u1, v)*dx - dt * dot(grad(p1), v)*dx
# Create variational formulation for k-eps
FK = (1. / dt) * dot(k - k0, phi)*dx \
+ dot(dot(u1, nabla_grad(k)), phi)*dx \
+ (nu + nut / 1.0) * inner(grad(k), grad(phi))*dx \
- dot(prod_k, phi)*dx + dot(react_k * k, phi)*dx
a_k = lhs(FK); l_k = rhs(FK)
FE = (1. / dt) * dot(e - e0, psi)*dx \
+ dot(dot(u1, nabla_grad(e)), psi)*dx \
+ (nu + nut / 1.3) * inner(grad(e), grad(psi))*dx \
- dot(prod_e, psi)*dx + dot(react_e * e, psi)*dx
a_e = lhs(FE); l_e = rhs(FE)
# main loop
iter_max = 1000
tol = 1e-5
for iter in range(iter_max):
# Solve N-S
A1 = assemble(a_1); b1 = assemble(l_1)
[bc.apply(A1, b1) for bc in bcu]
solve(A1, u1.vector(), b1)
A2 = assemble(a_2); b2 = assemble(l_2)
[bc.apply(A2, b2) for bc in bcp]
solve(A2, p1.vector(), b2)
A3 = assemble(a_3); b3 = assemble(l_3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u1.vector(), b3)
# Solve k-eps
AK = assemble(a_k); bK = assemble(l_k)
[bc.apply(AK, bK) for bc in bck]
solve(AK, k1.vector(), bK)
AE = assemble(a_e); bE = assemble(l_e)
[bc.apply(AE, bE) for bc in bce]
solve(AE, e1.vector(), bE)
# Make sure k and eps are positive
low_bound = 1e-16
dimension = len(k1.vector().get_local())
k1.vector()[:] = np.max([k1.vector()[:], low_bound * np.ones(dimension)], axis=0)
e1.vector()[:] = np.max([e1.vector()[:], low_bound * np.ones(dimension)], axis=0)
# Check if converged
err = u0 - u1
error = assemble(err**2*dx(mesh))
if error < tol:
print('Fully developed flow has been reached in %g iterations' % iter)
break
else:
print('iteration: %g, L2 error = %.6f (%.6f required)' % (iter+1, error, tol))
# Update all functions
u0.assign(u1)
k0.assign(k1)
e0.assign(e1)
plotting_enabled = 1
if plotting_enabled == 1:
# Plot solutions
cu = plot(sqrt(dot(u1, u1)), title = 'velocity')
plt.colorbar(cu)
plt.show()
cp = plot(p1, title = 'pressure')
plt.colorbar(cp)
plt.show()
ck = plot(k1, title = 'turbulent kinetic energy (k)')
plt.colorbar(ck)
plt.show()
ce = plot(e1, title = 'epsilon')
plt.colorbar(ce)
plt.show()
Thank you for reading,
I will highly appreciate anyone who takes their time and might lend me a hand with this 
.