Dear all,
Firstly I would like to make clear that I don't want to use this community to debug code! I just would like to clear some doubts. However, if my question is not well posed, or should not be asked here please, feel free to point this out and I will politely delete this post.
This is a follow up of this question which was already kindly answered.
https://fenicsproject.org/qa/13759/coupling-scalar-fields-yielding-dirichilet-expecting-function
I am trying to solve the following system of nonlinear time-dependent PDE's:
Those problems were solved and I was finally able to run the following code:
from fenics import *
import numpy as np
T = 100.0 # final time
num_steps = 1000 # number of time steps
dt = T / num_steps # time step size
L = 0.12
# Create mesh and define function space
nx = 20
mesh = IntervalMesh(nx,0,L)
P1 = FiniteElement('P', interval , 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)
# Definition of physical properties
k = 1.67
rho = 2200
c_p = 1100
B_T = Constant('1')
B_P = Constant('0.000001')
C_a = lambda T: 3.5 * 1e5 * (374.15 - T)**(1 / 3)
C_w = 4100
# Definition of coefficients (They will depend on T and P on further codes)
def A1(T,P):
return 0.15
def A2(T,P):
return - C_a(T) * 0.15
def A3(T,P):
return -0.3
def A4(T,P):
return rho * c_p + C_a(T) * 0.3
def A5(T,P):
return 1e-13
def A6(T,P):
return - C_w * 1e-13
# Define boundary conditions
# Dirichilet
T_D = Expression('t * 1 + 20.0', degree=2, t=0)
bcl_T = 'near(x[0],0)'
bc = DirichletBC(V.sub(0), T_D, bcl_T)
# Robin BC
class bc_right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0.12)
class bc_left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0.0)
bcr_P = bc_right()
bcr_T = bc_right()
bcl_P = bc_left()
boundaries = FacetFunction('size_t', mesh)
boundaries.set_all(0)
bcr_P.mark(boundaries,1)
bcl_P.mark(boundaries,1)
bcr_T.mark(boundaries,2)
ds = ds(subdomain_data=boundaries)
# Define initial values
T_0 = Constant('20')
P_0 = Constant('2900')
T_n = interpolate(T_0, V.sub(0).collapse())
P_n = interpolate(P_0, V.sub(1).collapse())
# Define variational problem
u = Function(V)
v_2, v_1 = TestFunctions(V)
T, P = split(u)
T_en = Constant('20') #Ambient temperature for robin condition
P_en = Constant('2090') #Ambient pressure for robin condition
f = Constant('0') #Zero source
F = (A1(T,P)*P+A3(T,P)*T)*v_1*dx - (A1(T,P)*P_n+A3(T,P)*T_n)*v_1*dx - dt * A5(T,P)* dot(grad(P), grad(v_1))*dx + \
dt * B_P * (P - P_en) * v_1 * ds(1) + (A2(T,P)*P + A4(T,P)*T)*v_2*dx - (A2(T,P)*P_n + \
A4(T,P)*T_n)*v_2*dx - dt * k * dot(grad(T), grad(v_2))*dx -dt * A6(T,P) * inner(grad(T), grad(P)) * v_2 * dx -\
dt * B_T * (- T + T_en) * v_2 * ds(2) + dt * f * v_2 * dx
# Time-stepping
t = 0
values_vert_T = [] # lists to save the temperature and pressure values
values_vert_P = []
progress = Progress('Time-stepping')
set_log_level(PROGRESS)
for n in range(100):
# Update current time
t += dt
T_D.t = t
# Compute solution
J = derivative(F, u) #Jacobian
solve(F == 0, u, bcs=bc, J=J)
T_, P_ = u.split(deepcopy=True)
# Plot solution
values_vert_T.append(T_)
values_vert_P.append(P_)
# Update previous solution
T_n.assign(T_)
P_n.assign(P_)
progress.update(t / num_steps)
However, after running this code, my solution presents big oscilations both in Temperature and in Pressure.
Would this be a hint that I made a mistake on the variational formulation? Or I should try to use different solvers?
Thank you for the help!
All the best, Murilo.
