### Tabulated material properties in multiple domains

45
views
0
9 weeks ago by
Hello,

I am doing simulations of heat and moisture transfer in building wall. The problem is governed by relatively simple diffusion equations but the material properties have a complex relation to the solution and are often given in tabular form. The wall typically consists of multiple materials.

Below you will see a test program for one-dimensional heat conduction with only one material. I have defined a Function lambdat for heat conductance and use
def build_lambdat_interpolator():
""" Use linear interpolation for tabulated lambdat values.

A constant value is used for extrapolation.
lambdat
|  ____
|      \
|       \_____  <- constant for extrapolation
|____________________
T
"""

# This table would be red from a file in real implementatation
table = sp.array([[0,0.01],
[5,0.05],
[10,0.10],
[15,0.15],
[20,0.20],
])
table[:,0] += C2K
return interp1d(table[:,0],table[:,1],
bounds_error=False,
fill_value=(table[0,1],table[-1,1])
)
lambdat_interpolator = build_lambdat_interpolator()

def update_lambdat(T, lambdat):
"""Update thermal conductivity
"""
Tar = T.vector().get_local() # .get_local() not nesessary
lambdat.vector().set_local(lambdat_interpolator(Tar))​

to set new heat conductivity values using interpolation from Scipy. It works just fine but extending this method to 3D with multiple materials seems tedious.

What is the best practice for this type of problem?

Full code below
from __future__ import division, print_function
import fenics as fe
from matplotlib import pyplot as plt
import scipy as sp
from scipy.interpolate import interp1d

C2K = 273.15

# Parameters
L       = 0.3
n       = 10
T_in    = 20 + C2K
T_out   = 0 + C2K

# Picard iteration
maxiter = 100
tol     = 1e-9

#############################################################################
# PRE
#############################################################################

# Create mesh and define function space
mesh = fe.IntervalMesh(n, 0, L)

# NOTE! In the code it is assumed that we are using linear elements!
# Something might break with higher order elements. Or maybe not.
# It is simply not tested or thought trhough.
V = fe.FunctionSpace(mesh, "CG", 1)

## Boundaries
inside_boundary_left   = fe.CompiledSubDomain('on_boundary && near(x[0], 0)')
outside_boundary_right = fe.CompiledSubDomain('on_boundary && near(x[0], L)',
L=L)

# Fixed temperature inside
T_in_const = fe.Constant(T_in)
bc_in  = fe.DirichletBC(V, T_in_const, inside_boundary_left)
# Fixed temperature outside
T_out_const = fe.Constant(T_out)
bc_out = fe.DirichletBC(V, T_out_const, outside_boundary_right)
bc = [bc_in, bc_out]

#############################################################################
# Material properties
#############################################################################
def build_lambdat_interpolator():
""" Use linear interpolation for tabulated lambdat values.

A constant value is used for extrapolation.
lambdat
|  ____
|      \
|       \_____  <- constant for extrapolation
|____________________
T
"""

# This table would be red from a file in real implementatation
table = sp.array([[0,0.01],
[5,0.05],
[10,0.10],
[15,0.15],
[20,0.20],
])
table[:,0] += C2K
return interp1d(table[:,0],table[:,1],
bounds_error=False,
fill_value=(table[0,1],table[-1,1])
)
lambdat_interpolator = build_lambdat_interpolator()

def update_lambdat(T, lambdat):
"""Update thermal conductivity
"""
Tar = T.vector().get_local() # .get_local() not nesessary
lambdat.vector().set_local(lambdat_interpolator(Tar))

#############################################################################
# SOLVER
#############################################################################

# HEAT SOLVER
T = fe.TrialFunction(V)
v = fe.TestFunction(V)
T_old = fe.interpolate(T_in_const, V)

# Heat conductivity
lambdat = fe.interpolate(T_in_const, V)
update_lambdat(T_old, lambdat)

# Weak form
a, L = fe.lhs(F), fe.rhs(F)
T = fe.Function(V)

# Picard iteration
eps = sp.inf
it = 0
while eps > tol and it < maxiter:
it += 1
fe.solve(a == L, T, bc)
diff = T.vector().get_local() - T_old.vector().get_local()
eps = sp.linalg.norm(diff, ord=sp.Inf)
print("it", it, "eps",eps)

# Update material properties
update_lambdat(T, lambdat)
T_old.assign(T)

#############################################################################
# POST
#############################################################################

## Plot solution and mesh
plt.figure(0)
fe.plot(mesh, title="Mesh")
plt.figure(1)
fe.plot(T-C2K, linestyle="-", color="r",title="T")

​
Community: FEniCS Project