Singularities in heat conduction with inhomogeneous properties

General
1

Dear all,

I am working on the phase field model to generate cracks due to thermal shock. The basic idea is throwing a hot plate into cool water, then we would see the cracks on the hot plate. In the phase model, the thermal conductivity would decrease at the location of cracks. The equation is k = (1 - damage)^2 * k0. damage stands for the crack, when damage = 1, that location is completely a crack.

The problem is that, at the cracks, with the heat conduction keeps going, I will have a few points with temperature far larger than the initial temperature (680 K) gradually. Please see the pictures below. Picture1 is the temperature field of the whole domain. Picture 2 is the crack pattern. Picture 3 is the singularities in one crack branch.

pic1
pic11172×324 40.1 KB

pic2
pic21174×331 57.4 KB

pic3
pic31181×644 41.6 KB

I made a simple code of the crack problem. Initially, I will set a part of the domain as crack with damage = 1 and with a fixed width = n * meshsize. Then, I will do heat conduction in the domain with initial temperature as 680 K, top and right boundaries as insulation, and bottom and left boundaries as 300 K. I found that, if the crack width is small, the maximum temperature (generally at the bottom of the crack) in the domain would be larger than 680 K, and then gradually decrease. If the crack width is large, the maximum temperature would keep going to a very large value. The code is attached as below. Can any one help me with this problem?

import numpy as np
from mpi4py import MPI
import ufl
from ufl import ds, dx, grad, inner, VectorElement, FacetNormal, SpatialCoordinate, TrialFunction, TestFunction, Measure
from dolfinx import mesh, fem, io, plot, default_scalar_type
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc
import pyvista
from petsc4py import PETSc

initial_dt = 5e-9    
conductiontime = 5e-5
t = 0
tol = 1e-7
small = 1e-8

L = 0.025
H = 0.005
nx = 500
ny = 100
meshsize = L / nx

timestep = initial_dt

mesh_comm = MPI.COMM_WORLD
msh = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([L, H])], [nx, ny], mesh.CellType.triangle)

V = fem.FunctionSpace(msh, ("Lagrange", 1))
WW = fem.FunctionSpace(msh, ("DG", 0))
x = SpatialCoordinate(msh)
T = TrialFunction(V)
v = TestFunction(V)

# set the initial crack part with width = n * meshsize, in the crack, damage = 1, other wise, damage = 0
damage = fem.Function(V)
def damage_interpolation(x):
    if x[0] <= 0.01:
        return 0
    elif x[0] > (0.01 + 4 * meshsize): # set the width of the crack
                                       # if width >= 3 * meshsize, the max temperature will keep increasing
                                       # if width <= 2 * meshsize, the max temperature will be only larger than 680 K at the beginning steps
        return 0
    else:
        return 1

dof_coords = V.tabulate_dof_coordinates()
damage_values = np.array([damage_interpolation(x) for x in dof_coords], dtype=np.float64)
damage.vector.setArray(damage_values)
damage.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)

alpha = fem.Function(WW)

def alpha_PF(p):
    return ((1.0 - p)**2 + small) * 300 / 2450 /0.775

alpha_expr = fem.Expression(alpha_PF(damage), WW.element.interpolation_points())
alpha.interpolate(alpha_expr)
alpha.x.scatter_forward()

# alpha = 300 / 2450 /0.775  # if alpha is constant, everything will be all good

def temp_bc(x):
    return np.logical_or(np.isclose(x[0], 0), np.isclose(x[1], 0))

dofs_crack = fem.locate_dofs_geometrical(V, temp_bc)
temp_boundary = fem.Function(V)
temp_boundary.interpolate(lambda x: 300.0 + 0 * x[0])
temp_boundary.x.scatter_forward()
bcs = [fem.dirichletbc(temp_boundary, dofs_crack)]

dt = fem.Constant(msh, default_scalar_type(initial_dt))

T_1 = fem.Function(V)
T_1.interpolate(lambda x: 680.0 + 0 * x[0])
T_1.x.scatter_forward()
Th = fem.Function(V)

a =  inner(T, v) * dx + dt * inner(alpha * grad(T), grad(v)) * dx 
L = inner(T_1, v) * dx

bilinear_form = fem.form(a)
linear_form = fem.form(L)
A = fem.petsc.create_matrix(bilinear_form)
b = fem.petsc.create_vector(linear_form)

solver = PETSc.KSP().create(msh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.GMRES)
solver.getPC().setType(PETSc.PC.Type.JACOBI)
solver.setTolerances(rtol=1e-8)
solver.setTolerances(atol=1e-12)
solver.setTolerances(max_it=1000)


while (t<conductiontime):
    
    t += timestep
    
    A.zeroEntries()
    fem.petsc.assemble_matrix(A, bilinear_form, bcs=bcs)
    A.assemble()
    with b.localForm() as loc_b:
        loc_b.set(0)
    fem.petsc.assemble_vector(b, linear_form)
    fem.petsc.apply_lifting(b, [bilinear_form], [bcs])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    fem.petsc.set_bc(b, bcs)

    solver.solve(b, Th.vector)
    Th.x.scatter_forward()  

    T_1.x.array[:] = Th.x.array
    T_1.x.scatter_forward()    

    max_temp = np.max((Th.x.array))
        
    print('time step: ', round(t, 9), "max temp:", round(max_temp, 3))

    if  round(t, 9) * 1e8 % 1 == 0:   
        topo, types, geom = plot.vtk_mesh(V)
        grid = pyvista.UnstructuredGrid(topo, types, geom)
        grid.point_data["temperature"] = Th.x.array.real
        grid.set_active_scalars("temperature")
        plotter = pyvista.Plotter()
        plotter.add_text("temperature", position="upper_edge", font_size=14, color="black")
        plotter.add_mesh(grid, show_edges=False)
        plotter.view_xy()
        plotter.show()
2

After some investigations, I think a way to solve this problem is to put an insulation boundary condition between the crack and other part of the material. However, I don’t know how to do that.

Or some others may have a better way to solve this problem?

Thanks!

3

I don’t recognize the heat equation by looking at your weak form. Are you sure it is equivalent to the common heat equation if you have no crack?

Also, if k is discontinuous across space, wouldn’t this means that the PDE is non linear, in which case you should use a nonlinear solver instead of the one you’re using?

About your last question, insulation means no heat flux, it means that the temperature across the spatial direction is constant. So a Dirichlet b.c. setting a fixed temperature at the crack surface and corresponding material surface should work, I believe.

4

Hello raboynics,

Thanks for your reply.

  1. To my understanding, the heat equation should be the same with that has no cracks. The only thing we need to change is the thermal conductivity. But I might be wrong.
  2. This problem is like applying two materials in the domain. It should be similar to the tutorial case Defining subdomains for different materials, so I think it’s still a linear problem.
  3. Yes, by setting a fixed temperature at the crack surface should work, though this way could be a little complex. However, I don’t know how to do that.

Thanks a lot for your reply!

5

I have ran your code as is (with some modifs just for the visuals, I cannot seem to be able to use pyvista with a docker container without jupyter notebook). In case this helps somebody to debug your code, the max temperature seems to go higher and higher, even higher than 900 K after 6 microseconds.

crack
crack1920×1440 71.8 KB

Edit: If I change the degree of V to 2: V = fem.FunctionSpace(msh, ("Lagrange", 2)) and small to 1e-7 (so 10x bigger than yours), I get this:

crack
crack1920×1440 68.2 KB

6

Yes, correct. The maximum temperature would go larger and larger than the initial temperature (680 K). BTW, the maximum temperature is always located at the crack bottom, where k = 0. This means by setting inhomogeneous properties may introduce some problems at the boundary.

By decreasing the crack width, like by setting the width = 1 * meshsize, the maximum temperature would only be larger than 680 K at first a few microseconds, and then go back to 680 K very quickly.

7

I could be mistaken but I think this is incorrect. You’re right that a thermally-insulated boundary condition (i.e. no heat flux) has a constant temperature in the spatial direction (i.e. the derivative of temperature w.r.t. the boundary’s normal vector is zero), but that’s not the same as setting a fixed temperature. Instead you should apply a homogeneous Neumann condition at the boundary. Whereas Dirichlet boundary conditions can act as thermal energy sources/sinks, Neumann conditions do not behave like that and do not directly affect the total amount of thermal energy over the domain.

8

I also get the same behavior (goes towards 680 K from above quickly) if I set “small” to 1e-2.
The thermal conductivity of air is like 0.025. If you go several orders of magnitude below this, it’s as if you had a vacuum. By looking at alpha: ((1.0 - p) ** 2 + small) * 300 / 2450 / 0.775, I see 2 values.

1: where p equals 1 (crack), so alpha equals 1.6e-9 (i.e. probably high vacuum).
2: where p = 0 (everywhere else but cracks), alpha = 0.16 which is essentially glass as you pointed out.

Why do you have a squared exponent in (1-p)**2, it has no effect whatsoever?

Setting small = 1e-2 means alpha can be worth one tenth to one-fifteenth of the thermal conductivity of air, i.e. some vacuum. Which seems reasonable to me, and doesn’t yield that crazy results, I think.

9

I didn’t express myself well enough. By fixed temperature I did not mean to fix it through time. I meant to set it up manually equal at the boundaries crack/glass domain. This value would have to be modified at each time step. I also do not see how to do it properly. But mathematically, I think this is correct.

10

(1-p)**2 is an equation in the phase field model, because in phase field model, the crack is assumed to be continuous from p=1 to p=0. The p actually will change with the temperature and strain in the material. As this is only a simple version of the code, I only set a crack in the middle directly (p=1 or p=0).

11

Yes, correct. An insulation boundary would be better. Do you know how to apply that in the middle of the material?