Steady k-epsilon turbulence implementation

Question

Hello there,

I am getting an error when I try to solve the k-epsilon equations for fluid flow.

Whilst solving nonlinear variational problem, the following error comes up:

[0]PETSC ERROR: Floating point exception!
[0]PETSC ERROR: Infinite or not-a-number generated in norm!

I can't figure out where this is coming from. ANY help would be hugely appreciated!

I haven't yet figured out how to implement wall functions so to start with I am just imposing Dirichlet BCS.

Apart from that it is all pretty standard.

I attach my code below:

from dolfin import *

parameters["form_compiler"]["cpp_optimize"] = True
parameters["allow_extrapolation"] = True
parameters["num_threads"] = 8


rho = Constant(1.0)
nu = Constant(1.0)
Cmu = Constant(1.0)
C1 = Constant(1.0)
C2 = Constant(1.0)
sigma = Constant(1.0)
cb = Constant(1.0)
l0 = Constant(1.0)

# define mesh
mesh = UnitSquareMesh(100, 100)


# define functions
V = FunctionSpace(mesh, 'CG', 1)
W = MixedFunctionSpace([V,V,V,V,V])
w = Function(W)
(u, v, p, k, e) = split(w)
(vu, vv, vp, vk, ve) = TestFunctions(W)

# influx conditions
u_inf = Constant(1.0)
k_inf = Constant(cb*(u_inf **2))
e_inf = Constant(Cmu*(k_inf **(3/2))/l0)


# specify subdomains, declared in C++ for potential speed increase
ns = CompiledSubDomain("on_boundary && (x[1] < DOLFIN_EPS || x[1] > 1.0 - DOLFIN_EPS)")
out = CompiledSubDomain("x[0] > 1.0 - DOLFIN_EPS")
inf = CompiledSubDomain("x[0] < DOLFIN_EPS")

# define boundary conditions
# no slip
bc_nsu  = DirichletBC(W.sub(0), Constant(0.0), ns)
bc_nsv  = DirichletBC(W.sub(1), Constant(0.0), ns)
bc_nsk  = DirichletBC(W.sub(3), Constant(0.0), ns)
bc_nse  = DirichletBC(W.sub(4), Constant(0.0), ns)

# influx
bc_infu = DirichletBC(W.sub(0), u_inf, inf)
bc_infv = DirichletBC(W.sub(1), Constant(0.0), inf)
bc_infk = DirichletBC(W.sub(3), k_inf, inf)
bc_infe = DirichletBC(W.sub(4), e_inf, inf)
bc_inp = DirichletBC(W.sub(2), Constant(0.0), inf)

bcs = [bc_nsu, bc_nsv, bc_nsk, bc_nse, bc_infu, bc_infv, bc_infk, bc_infe, bc_inp]

# define variational form
nu_t = Cmu * (k**2) / e

Fu = (u*Dx(u,0)*vu + v*Dx(u,1)*vu - p*Dx(vu,0)/rho + (nu + nu_t)*inner(nabla_grad(u),nabla_grad(vu)) \
      - 2*Dx(nu_t,0)*Dx(u,0)*vu - Dx(nu_t,1)*(Dx(v,0)+Dx(u,1))*vu )*dx

Fv = (u*Dx(v,0)*vv + v*Dx(v,1)*vv - p*Dx(vv,1)/rho + (nu + nu_t)*inner(nabla_grad(v),nabla_grad(vv))\
      - 2*Dx(nu_t,1)*Dx(v,1)*vv - Dx(nu_t,0)*(Dx(v,0)+Dx(u,1))*vv )*dx

Fm = (Dx(u,0)*vp + Dx(v,1)*vp)*dx

Fk = (Dx(u*k,0)*vk + Dx(v*k,1)*vk + nu_t*Dx(k,0)*Dx(vk,0) + nu_t*Dx(k,1)*Dx(vk,1) \
      - nu_t*(2*Dx(u,0) **2 + (Dx(v,0)+Dx(u,1)) **2 + 2*Dx(v,1) **2)*vk + e*vk)*dx

Fe = (Dx(u*e,0)*ve + Dx(v*e,1)*ve + nu_t*Dx(e,0)*Dx(ve,0)/sigma + nu_t*Dx(e,1)*Dx(ve,1)/sigma \
      - nu_t * C1 * (2* Dx(u,0) **2 + (Dx(v,0)+Dx(u,1)) **2 + 2 * Dx(v,1) **2) * ve * e / k + 2*C2*(e **2) * ve / k)*dx

F = Fu+Fv+Fm+Fk+Fe


# find form of jacobian
J = derivative(F, w)


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

prm = solver.parameters['newton_solver']
prm['maximum_iterations'] = 100
prm['absolute_tolerance'] = 1E-3
prm['relative_tolerance'] = 1E-3

solver.solve()

graham.benham · Answer 1 · Dec 8, 2015

I managed to fix it myself by providing an initial guess and placing a conditional expression on the eddy viscosity nu_t so that it doesn't blow up when e = 0.

This still is only the trivial solution (dirichlet bcd everywhere!). I would love it if someone showed me how to implement wall functions!!!!

CODE

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

parameters["form_compiler"]["cpp_optimize"] = True
parameters["allow_extrapolation"] = True
parameters["num_threads"] = 8

rho = Constant(1.0)
nu = Constant(1.0)
Cmu = Constant(1.0)
C1 = Constant(1.0)
C2 = Constant(1.0)
sigma = Constant(1.0)
cb = Constant(1.0)
l0 = Constant(1.0)

define mesh

mesh = UnitSquareMesh(100, 100)

define functions

V = FunctionSpace(mesh, 'CG', 1)
W = MixedFunctionSpace([V,V,V,V,V])
w = Function(W)
w=interpolate(Expression(("x[0]+x[1]","x[0]+x[1]","x[0]+x[1]","x[0]+x[1]","x[0]+x[1]")),W)

(u, v, p, k, e) = split(w)
(vu, vv, vp, vk, ve) = TestFunctions(W)

bcs=DirichletBC(W,Constant((0.0, 0.0, 0.0, 0.0, 0.0)),DomainBoundary())

define variational form

tol=Constant(0.01)
nu_t = conditional( lt(abs(e),tol), project(Constant(0.0),V), Cmu*(k**2)/e)

Fu = (uDx(u,0)vu + vDx(u,1)vu - pDx(vu,0)/rho + (nu + nu_t)inner(nabla_grad(u),nabla_grad(vu)) \
- 2Dx(nu_t,0)Dx(u,0)vu - Dx(nu_t,1)(Dx(v,0)+Dx(u,1))vu )dx

Fv = (uDx(v,0)vv + vDx(v,1)vv - pDx(vv,1)/rho + (nu + nu_t)inner(nabla_grad(v),nabla_grad(vv))\
- 2Dx(nu_t,1)Dx(v,1)vv - Dx(nu_t,0)(Dx(v,0)+Dx(u,1))vv )dx

Fm = (Dx(u,0)vp + Dx(v,1)vp)*dx

Fk = ((Dx(uk,0)vk + Dx(vk,1)vk + nu_tDx(k,0)Dx(vk,0) + nu_tDx(k,1)Dx(vk,1) \
- nu_t(2(Dx(u,0)2)+(Dx(v,0)+Dx(u,1))2+2*Dx(v,1)**2)*vk + e*vk))*dx

Fk = ((Dx(uk,0)vk + Dx(vk,1)vk + nu_tDx(k,0)Dx(vk,0)+ nu_tDx(k,1)Dx(vk,1) -nu_t*Dx(u,1)**2*vk + e*vk))*dx

Fe = ((Dx(ue,0)ve + Dx(ve,1)ve + nu_tDx(e,0)Dx(ve,0)/sigma + nu_tDx(e,1)Dx(ve,1)/sigma \
-nu_tC1(2Dx(u,0)2+(Dx(v,0)+Dx(u,1))2+2Dx(v,1)2)vee/k + 2C2(e2)ve/k))dx

F = Fu+Fv+Fm+Fk+Fe

find form of jacobian

J = derivative(F, w)

solve

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

prm = solver.parameters['newton_solver']
prm['maximum_iterations'] = 100
prm['absolute_tolerance'] = 1E-3
prm['relative_tolerance'] = 1E-3

solver.solve()

Hi Graham,
May I ask have you found a way to implement wall functions?
Thanks a lot! — Hrunx, Feb 29, 2016