Hello everyone,
I am trying to implement a conjugate heat transfer problem. My problem domain looks like this,
Earlier, I have done this in a segregated way, which I solve the Navier-Stroke Equation using splitting method (as explained in the FEniCS tutorial) in the fluid domain for pressure and velocity, then I apply the solved velocity to the energy equation to solve for the temperature distribution over the whole domain. According to the following code, this works fine for me.
# Define function space
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 2)
# Fluid and Solid Properties
muF = 0.89 # dynamic viscosity
rhoF = 1000 # density
cpF = 4220 # heat capacity
cdF = 0.6071 # heat conductivity
alpF = cdF/(rhoF*cpF) # heat diffusivity
rhoS = 2700 # density
cpS = 20 # heat capacity
cdS = 896 # heat conductivity
alpS = cdF/(rhoF*cpF) # heat diffusivity
#Velocity BCs
bcu = [DirichletBC(V.sub(0), Constant(0), mf, 3),
DirichletBC(V.sub(1), Constant(0), mf, 3),
DirichletBC(V.sub(0), Constant(0), mf, 6),
DirichletBC(V.sub(1), Constant(0), mf, 6)]
#Pressure BCs
bcp = [DirichletBC(Q, Constant(8), mf, 4),
DirichletBC(Q, Constant(0), mf, 5)]
#Temperature BCs
bcT = [DirichletBC(Q, Constant(20), mf, 4)]
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
TE = TrialFunction(Q)
Z = TestFunction(Q)
# Define functions for solutions at previous and current time steps
u_n = Function(V) #this is u at time n
u_ = Function(V) #this is u at time n+1
p_n = Function(Q) #this is p at time n
p_ = Function(Q) #this is p at time n+1
T_n = Function(Q) #this is T at time n
T_ = Function(Q) #this is T at time n+1
# Setting the control parameters
T = 20 # final time
num_steps = 400 # number of time steps
dt = T / num_steps # time step size
# Define expressions used in variational forms
U = 0.5*(u_n + u)
n = FacetNormal(mesh)
f = Constant((0, 0)) #source term
k = Constant(dt)
muF = Constant(muF)
rhoF = Constant(rhoF)
alpF = Constant(alpF)
rhoS = Constant(rhoS)
cdS = Constant(cdS)
alpS = Constant(alpS)
qS = 475000/cdS
# set initial condition for temperature
Tint = Constant(25)
T_n = interpolate(Tint, Q)
# Define strain-rate tensor
def epsilonF(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigmaF(u, p):
return 2*muF*epsilonF(u) - p*Identity(len(u))
# Define variational problem for step 1
F1 = rhoF*dot((u - u_n) / k, v)*dx(1) \
+ rhoF*dot(dot(u_n, nabla_grad(u_n)), v)*dx(1) \
+ inner(sigmaF(U, p_n), epsilonF(v))*dx(1) \
+ dot(p_n*n, v)*ds(3) - dot(muF*nabla_grad(U)*n, v)*ds(3) \
+ dot(p_n*n, v)*ds(4) - dot(muF*nabla_grad(U)*n, v)*ds(4) \
+ dot(p_n*n, v)*ds(5) - dot(muF*nabla_grad(U)*n, v)*ds(5) \
+ dot(p_n*n, v)*ds(6) - dot(muF*nabla_grad(U)*n, v)*ds(6) \
- dot(f, v)*dx(1)
a1 = lhs(F1)
L1 = rhs(F1)
# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx(1)
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx(1) - (rhoF/k)*div(u_)*q*dx(1)
# Define variational problem for step 3
a3 = dot(u, v)*dx(1)
L3 = dot(u_, v)*dx(1) - (k/rhoF)*dot(nabla_grad(p_ - p_n), v)*dx(1)
# Define variational problem for step 4
a4 = TE*Z*dx + k*dot(u_, grad(TE))*Z*dx + k*alpF*dot(grad(TE), grad(Z))*dx
L4 = T_n*Z*dx + alpS*k*qS*Z*ds(9) #q4*v*ds(4)
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
A4 = assemble(a4)
# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
[bc.apply(A4) for bc in bcT]
# Time-stepping
t = 0
for n in range(num_steps):
# Update current time
t += dt
# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1)
# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2)
# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3)
# Step 4: Calculate Temperature Distribution
b4 = assemble(L4)
[bc.apply(b4) for bc in bcT]
solve(A4, T_.vector(), b4)
# Update previous solution
u_n.assign(u_)
p_n.assign(p_)
T_n.assign(T_)
However, for some reasons as well as the above code take very long time until the solution converge, I need to change the implementation to the monolithic way, which I combine the steady momentum equation, continuity equation and energy equation together. This works when I use only a single domain (e.g. flow problem with temperature boundary conditions), unfortunately, it does not work when I apply to the same domain as shown in the above picture. In this way, I use nonlinear solver and I always encounter the diverge residual. The following code is the code that I encounter the problem. Can anyone suggest on what is wrong with my code?
# Define function spaces
V = VectorElement('P', triangle, 2)
P = FiniteElement('P', triangle, 1)
Q = FiniteElement('P', triangle, 1)
W = FunctionSpace(mesh, MixedElement([V, P, Q]))
# Define variational problem
(u, p, T) = TrialFunctions(W)
(v, q, Z) = TestFunctions(W)
bcs = [DirichletBC(W.sub(0).sub(0), Constant(0), mf, 3),
DirichletBC(W.sub(0).sub(1), Constant(0), mf, 3),
DirichletBC(W.sub(0).sub(0), Constant(0), mf, 6),
DirichletBC(W.sub(0).sub(1), Constant(0), mf, 6),
DirichletBC(W.sub(1), Constant(80), mf, 4),
DirichletBC(W.sub(1), Constant(0), mf, 5)]
# Surface normal
n = FacetNormal(mesh)
# Define parameters used in the weak form
mu = 0.89
rho = 1000
k = 0.6071
cp = 4220
alp = k/rho/cp
muS = 1
rhoS = 3000
kS = 1.5
cpS = 220
alpS = kS/rhoS/cpS
Fl = 0.00005
# Define initial conditions
e_u0 = Constant((0,0))
e_p0 = Constant(0)
e_T0 = Constant(0)
u0 = interpolate(e_u0, W.sub(0).collapse())
p0 = interpolate(e_p0, W.sub(1).collapse())
T0 = interpolate(e_T0, W.sub(2).collapse())
# Solution vectors
w = Function(W)
assign(w, [u0, p0, T0])
(u, p, T) = split(w)
# Weak form of the steady state fluid problem
F = rho*inner(grad(u)*u, v)*dx(1) \
+ mu*inner(grad(u), grad(v))*dx(1) \
+ inner(grad(p),v)*dx(1) \
- div(u)*q*dx(1) \
+ inner(dot(u,grad(T)),Z)*dx(1) \
+ alp*inner(grad(T), grad(Z))*dx(1) + alpS*inner(grad(T), grad(Z))*dx(2) \
- 1000*Fl*Z*ds(9) \
- mu*dot(grad(u)*n,v)*ds(3) \
- mu*dot(grad(u)*n,v)*ds(4) \
- mu*dot(grad(u)*n,v)*ds(5) \
- mu*dot(grad(u)*n,v)*ds(6)
# Jacobian matrix to be supplied to the non linear solver
J = derivative(F, w) # Jacobian
# Call the non'linear solver
solve(F == 0, w, bcs, J=J)
Thank you very much.