I am trying to solve the heat flow equation near the melting point. As a consequence, if the temperature exceeds the melting point, the material becomes liquid and the thermal conductivity is different. Is it possible to create a subdomain for the liquid and other for the solid? These subdomains would change after each time step depending on the temperature at each grid point.
Thanks in advance!
Assuming you’re solving the same equation with a different conductivity in the two domains (and not, e.g., changing the partial differential equation system by introducing convection in the liquid region), the most expedient approach would probably be to define the thermal conductivity using a UFL conditional expression, as follows:
k = conditional(gt(T,T_melt),k_liquid,k_solid)
where T is the unknown temperature, T_melt is the melting point, k_liquid is the thermal conductivity when T is greater than (gt) the melting point, and k_solid is the conductivity when T is below (or equal to) the melting point. (The two coefficients k_liquid and k_solid could also be themselves functions of T defined earlier.)
Assuming T is the current/unknown temperature, the above would introduce the change in conductivity implicitly, and require a nonlinear solve at each time step. You could also try replacing T in the conditional with the previous time step’s solution to determine the subdomains explicitly, which would be faster and keep the problem linear (assuming k_liquid and k_solid are independent of T or also defined explicitly using the previous time step’s solution), but this may be less stable at large time steps.
I tried this:
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
#kappa = 31.3
kappa = conditional(gt(u,Constant(1500)),31,11)
But I get this error:
Found Argument in , this is an invalid expression.
Traceback (most recent call last):
File “/scratch/snx3000/jfernnde/fenics/laser_tprofile/run_cscs.py”, line 177, in
SolveProblem(mesh,boundaries,subdomains,t)
File “/scratch/snx3000/jfernnde/fenics/laser_tprofile/run_cscs.py”, line 153, in SolveProblem
a, R = lhs(F), rhs(F)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/formoperators.py”, line 79, in lhs
return compute_form_lhs(form)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/formtransformations.py”, line 379, in compute_form_lhs
return compute_form_with_arity(form, 2)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/formtransformations.py”, line 345, in compute_form_with_arity
return map_integrands(_transform, form)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/map_integrands.py”, line 38, in map_integrands
mapped_integrals = [map_integrands(function, itg, only_integral_type)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/map_integrands.py”, line 38, in
mapped_integrals = [map_integrands(function, itg, only_integral_type)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/map_integrands.py”, line 46, in map_integrands
return itg.reconstruct(function(itg.integrand()))
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/formtransformations.py”, line 341, in _transform
e, provides = pe.visit(e)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/transformer.py”, line 110, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/transformer.py”, line 110, in
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/transformer.py”, line 114, in visit
r = h(o)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/formtransformations.py”, line 275, in linear_indexed_type
part, provides = self.visit(expression)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/transformer.py”, line 110, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/transformer.py”, line 110, in
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/transformer.py”, line 114, in visit
r = h(o)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/formtransformations.py”, line 275, in linear_indexed_type
part, provides = self.visit(expression)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/transformer.py”, line 114, in visit
r = h(o)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/formtransformations.py”, line 275, in linear_indexed_type
part, provides = self.visit(expression)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/transformer.py”, line 110, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/transformer.py”, line 110, in
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/transformer.py”, line 114, in visit
r = h(o)
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/algorithms/formtransformations.py”, line 73, in expr
error(“Found Argument in %s, this is an invalid expression.” % ufl_err_str(x))
File “/users/jfernnde/bin/miniconda3/envs/fenics/lib/python3.9/site-packages/ufl/log.py”, line 172, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Found Argument in , this is an invalid expression.
Do you know what’s wrong?
Once you have a piecewise function (i.e. your conditional) which is a function of the solution, your problem becomes nonlinear. I’d recommend consulting the nonlinear Poisson demo.
The error is telling you that a TrialFunction Argument has been found where it shouldn’t.
But that example is not given as a function of time, right? Can it be used in that case? I did not see complex models like this one