Anisotropic properties with different materials

1

I am simulating heat diffusion through several layers of dissimilar materials, with one of them having an anisotropic thermal conductivity.
I defined the thermal conductivity function as:

class Thermal_Conductivity(UserExpression):
    
    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs) # This part is new!
        self.markers = markers
    
    def eval_cell(self, values, x, cell):
        for mat in materials
            im = materials.index(mat)
            if self.markers[cell.index] == association_table[mat]:
                kx = kx_materials[im]
                ky = ky_materials[im]
                kz = kz_materials[im]
                values = as_tensor([[kx, 0., 0.],[0., ky, 0.],[0., 0., kz]])
                
k = Thermal_Conductivity(subdomains, degree=1)

The solution converges but yields singularities with very hot spots, even with kx=ky=kz.
Is this the correct way to do it?

2

Edit:
I think this post advising to use the curl instead of the gradient might be the solution: Anisotropic material definition and results issue - #2 by kamensky
However even with curl I still get NaNs in my solution…

3

Can you post a MWE? I’m not convinced your code is written correctly since you reassign values inside your eval_cell method, rather than populating its members.

Also you’re using as_tensor which is a ufl function, not a numerical operation.

And finally, you haven’t specified the shape of the tensor, you haven’t implemented the value_shape() method.

1 Like
4

Thank you for these hints. I have tried adding value_shape() but it yielded an error with the variational form.
Therefore I produced a MWE as suggested and uploaded it here (including the mesh): Fenics/anisotropy.zip at 0f7845fdf94b00e990af70632c401925d438243f · bragostin/Fenics · GitHub
The python code itself is:

import numpy as np
import matplotlib.pyplot as plt
from msh2xdmf import import_mesh, msh2xdmf
from dolfin import *
from ufl import Min

# LOAD MESH
mesh, boundaries, subdomains, association_table = import_mesh(prefix='anisotropy', dim=3, subdomains=True)

# Boundary conditions
HTC = 10000. # (W/m2 K)
Tref = 60. # (°C)
qc = 100. * 1e4 # (W/m2)

# Define subdomain markers
markers = MeshFunction('size_t', mesh, 3, mesh.domains())
print("association_table = ", association_table)

# Define thermal conductivity
class Thermal_Conductivity(UserExpression):
    
    # Returns error : ufl.log.UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (3, 3) and free indices with labels ().
    #def value_shape(self):
    #    return (3,3)
    
    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs) # This part is new!
        self.markers = markers
    
    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == association_table['base']:
            #values[0] = 100. # isotropic --> works
            values = as_tensor([[100., 0, 0],[0, 100., 0],[0, 0, 100.]])
        else:
            #values[0] = 10. # isotropic --> works
            values = as_tensor([[100., 0, 0],[0, 100., 0],[0, 0, 100.]])
            
k = Thermal_Conductivity(subdomains, degree=1)

# Build function space
V = FunctionSpace(mesh, 'CG', 1)
T = TrialFunction(V)
s = TestFunction(V)
n = FacetNormal(mesh)
ds_bc = ds(subdomain_data=boundaries)

bcs = []

# DEFINE VARIATIONAL FORM
F = (
    - k*inner(grad(T), grad(s))*dx
    - HTC*(T - Tref)*s*ds_bc(association_table["htc"])
    + qc*s*ds_bc(association_table["top1"])
    )

a, L = lhs(F), rhs(F)

# Create VTK files for visualization output
vtkfile_T = File('Temperature.pvd')

T = Function(V, name='Temperature')
#solve(a == L, T, bcs)
solve(a == L, T, bcs, solver_parameters={'linear_solver':'mumps'})

# Save solution to file (VTK)
vtkfile_T << T

# Isotropic max temperature = 310 --> OK
# Anisotropic max temperature = 12000 --> not OK
5

A couple of things indicating your formulation is incorrect:

  • k*inner(grad(T), grad(s))*dx doesn’t make sense for the heat equation’s variational form if you want k to be a tensor. This is why UFL was throwing an error when you set the shape to (3, 3)
  • values = as_tensor([[100., 0, 0],[0, 100., 0],[0, 0, 100.]]) You’re overwriting the values variable here which makes the underlying callback redundant. Check the demos.

See the following code which gives me something that (at least) looks reasonable.

import numpy as np
from dolfin import *

mesh = UnitSquareMesh(32, 32)


# Define thermal conductivity
class Thermal_Conductivity(UserExpression):

    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs) # This part is new!
        self.markers = markers
    
    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == 0:
            k = np.array([[1.0, 0.0],
            			  [0.0, 1.0]])
            values[:] = k.flatten()
        else:
            k = np.array([[2.0, -1.0],
            			  [-1.0, 2.0]])
            values[:] = k.flatten()

    def value_shape(self):
    	return (2, 2)
            

subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
AutoSubDomain(lambda x, on: x[0] > 0.5).mark(subdomains, 1)

k = Thermal_Conductivity(subdomains, degree=1)


# Build function space
V = FunctionSpace(mesh, 'CG', 1)
T = TrialFunction(V)
s = TestFunction(V)
n = FacetNormal(mesh)

bcs = [DirichletBC(V, Constant(1.0), "on_boundary")]

f = Constant(1.0)
F = inner(k*grad(T), grad(s))*dx - f*s*dx

a, L = lhs(F), rhs(F)

# Create VTK files for visualization output
vtkfile_T = File('Temperature.pvd')

T = Function(V, name='Temperature')
solve(a == L, T, bcs, solver_parameters={'linear_solver':'mumps'})

# Save solution to file (VTK)
vtkfile_T << T

image

2 Likes
6

Thank you very much! That works perfectly and yields physical solutions.
I noticed that with an anisotropic thermal conductivity, the solver takes about 4 times longer to converge. Is that to be expected?
Also, k = Thermal_Conductivity(subdomains, degree=0) instead of degre=1 yields the same solution but is faster.