Hello everyone. I am trying to model heat transfer with two domains but it seems that the inner domain has its temperature kept constant for some reason ( see picture) . It the computed example it should not be the case as the physical properties of the two domain are the same.
Here is the code for a MWE
from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace,
form, locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem, assemble_matrix, create_vector, assemble_vector, apply_lifting, set_bc
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import create_unit_square, locate_entities, locate_entities_boundary
from ufl import (SpatialCoordinate, TestFunction, TrialFunction,
dx, grad, inner, Measure, dot, jump, FacetNormal, avg, FiniteElement, MixedElement)
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import time
def project(function, space):
p = TrialFunction(space)
q = TestFunction(space)
a = inner(p, q) * dx
L = inner(function, q) * dx
problem = LinearProblem(a, L, bcs = [])
return problem.solve()
###############################################################################
################################# PARAMETER ###################################
###############################################################################
# PhysicalProperties
# density (kg.m^-3)
rho_1 = 8000.0e-6
rho_2 = 8000.0e-6
# capacity (J/Kg.K)
cp_1 = 435.0e6
cp_2 = 435.5e6
# conductivity (W/m.K)
k_1 = 10.0e3
k_2 = 10.0e3
# Boundary condtions ---------------------------------------------------------
# build plate temperature (K)
T_plate = 353.0
# initial temperature (K)
T0 = 293.0
T0_1 = 293.0
T0_2 = 320.0
# Time ----------------------------------------------------------------------
TotalTime = 1000.0
t0=time.time()
num_steps = 50 # number of time steps
dt = TotalTime / num_steps # time step size
t = 0 # current time
CYL_MARK = 2
EXT_MARK = 1
FACE_MARK = 1
###############################################################################
################################# Mesh ########################################
###############################################################################
with XDMFFile(MPI.COMM_WORLD, "CYL_IN_BOX_SLICE.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
with XDMFFile(MPI.COMM_WORLD, "facet_CYL_IN_BOX_SLICE.xdmf", "r") as xdmf:
ft = xdmf.read_meshtags(mesh, name="Grid")
ds = Measure('ds', domain=mesh, subdomain_data=ft)
###############################################################################
################################# MULTIDOMAIN #################################
###############################################################################
RhoCp_values = [rho_1*cp_1,rho_2*cp_2]
Conductivity_values = [k_1,k_2]
Tinit_values = [T0_1,T0_2]
Q = FunctionSpace(mesh, ("DG", 0))
Conductivity = Function(Q)
for i in range(len(Conductivity_values)):
cells = ct.find(i+1)
Conductivity.x.array[cells] = np.full_like(cells, Conductivity_values[i], dtype=default_scalar_type)
RhoCp = Function(Q)
for i in range(len(Conductivity_values)):
cells = ct.find(i+1)
RhoCp.x.array[cells] = np.full_like(cells, RhoCp_values[i], dtype=default_scalar_type)
T_INIT = Function(Q)
T_INIT.name = "Tinit"
for i in range(len(Tinit_values)):
cells = ct.find(i+1)
T_INIT.x.array[cells] = np.full_like(cells, Tinit_values[i], dtype=default_scalar_type)
###############################################################################
################################# FORMULATION #################################
###############################################################################
# Define variational problem
V = FunctionSpace(mesh, ('CG', 1))
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(mesh, PETSc.ScalarType(0))
n = FacetNormal(mesh)
# initialization
Tinit = Function(V)
Tinit.name = "Tinit"
# variable, need to be projected form Q onto V
Tinit = project(T_INIT,V)
Tinit.name = "Tinit"
xdmf = XDMFFile(mesh.comm, "./RES/Temperature.xdmf", "w")
xdmf.write_mesh(mesh)
Tinit = project(T_INIT,V)
Tinit.name = "Tinit"
Timp_facets = ft.find(FACE_MARK)
Timp_dofs = locate_dofs_topological(V, mesh.topology.dim - 1, Timp_facets)
bc = dirichletbc(default_scalar_type(T_plate), Timp_dofs, V)
a = RhoCp*u*v*dx
a += dt*dot(Conductivity*grad(u), grad(v))*dx
L = (RhoCp*Tinit + dt * f) * v * dx
bilinear_form = form(a)
linear_form = form(L)
A = assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = create_vector(linear_form)
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
# Define solution variable, and interpolate initial solution for visualization in Paraview
sol = Function(V)
sol.name = "Temperature"
sol.interpolate(Tinit)
for i in range(num_steps):
t += dt
print("current time is " + str(t) + " s")
print("Completion is " + str(100*(t/TotalTime))+ " %")
# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, linear_form)
# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [bilinear_form], [[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bc])
# Solve linear problem
solver.solve(b, sol.vector)
sol.x.scatter_forward()
# Update solution at previous time step (u_n)
Tinit.x.array[:] = sol.x.array
# Write solution to file
xdmf.write_function(sol, t)
xdmf.close()
Not sure if I am missing something …
Any advice would be very appreciated .
Thanks,
Quentin
