Defining heat flux (neumann) boundary condition on centroids of boundary panels

dolfinx
,
1

Hello!

Background info about simulation

I’ve been trying to setup a simulation which solves the following heat equation:

\frac{\partial u(\vec{X},t)}{\partial t} = \nabla \cdot \left( \alpha(\vec{X},u) \nabla u(\vec{X},t)\right) + f(\vec{X},t) \; \; \; \; \text{in} \; \Omega \\

with the initial condition

u(\vec{X},t_0) = u_0(\vec{X}) \; \; \; \; \; \text{in} \; \Omega \\

and the following Neumann boundary condition:

-k(\vec{X},u) \frac{\partial u}{\partial n}(\vec{X},t) = q(\vec{X},t) \; \; \; \; \text{on} \; \partial \Omega

where
u(\vec{X},t) is the temperature of the cell at location \vec{X} and time t [K].
\alpha(\vec{X},u) is the thermal diffusivity of the material at location \vec{X} and temperature u [m^2/s].
f(\vec{X},t) is the increase in temperature at location \vec{X} and time t due to internal processes [K].
k(\vec{X},u) is the thermal conductivity of the material at location \vec{X} and temperature u [W/mK]
q(\vec{X},t) is the heat flux at location \vec{X} and time t [W/m^2].

for this I’ve derived following variational formulation:

a(u,v) = \int_\Omega \left( u^{n+1} v - \Delta t \alpha \nabla u^{n+1} \cdot \nabla v \right) d\nu
L_{n+1}(v) = \int_\Omega \left( u^n v + \Delta t f^{n+1} v\right) d\nu + \int_{\partial \Omega} \Delta t \alpha \frac{- q}{k} v ds

Description of issue

I am struggling to understand how to define the function q(\vec{X},t) using dolfinx given the input data I have: From a trajectory propagation and using an engineering model (Quinn 2000) I obtain the heat flux on each boundary panel of an arbitrary mesh. In other words: I have a value for q at the centroids of each boundary panel. I have not found a way to pass this data into the simulation. I have looked at creating a function space which is only defined on the centroids of the boundary panels but have not fully understood how to do this (if this is even possible). I would then create a function in this function space to which I would write my values for q.

Here is a MWE, I’ve marked the places where I would need a function space which is only defined on the centroids of the boundary panels with # TODO: (line 50 and 89)

import numpy as np
import ufl
from dolfinx import fem, mesh, io, plot
import pyvista
from petsc4py import PETSc
from dolfinx.fem.petsc import assemble_vector, assemble_matrix, create_vector
import matplotlib as mpl
from mpi4py import MPI


cubeMesh = mesh.create_unit_cube(MPI.COMM_WORLD, 10, 10, 10)

def boundary_definition_func(X):
        
    out = np.sum(np.array([
        np.isclose(X[0], 0),
        np.isclose(X[0], 1),
        np.isclose(X[1], 0),
        np.isclose(X[1], 1),
        np.isclose(X[2], 0),
        np.isclose(X[2], 1),
    ]), axis=0) > 0
            
    return out 

boundary_definition_func.tag = 1

externalPanelRegionTags = [boundary_definition_func.tag]
externalPanelTagsList = [mesh.locate_entities(cubeMesh, cubeMesh.topology.dim-1, boundary_definition_func)]
    
externalPanelMarkers = []
for regionMarker, panelTags in zip(externalPanelRegionTags, externalPanelTagsList):
    externalPanelMarkers.append(np.full_like(panelTags, regionMarker))
    
allBoundaryPanelTags = np.hstack(externalPanelTagsList).astype(np.int32)
allBoundaryPanelRegionTags = np.hstack(externalPanelMarkers).astype(np.int32)

boundaryPanelMeshTagsObj = mesh.meshtags(cubeMesh, cubeMesh.topology.dim-1, allBoundaryPanelTags, allBoundaryPanelRegionTags)

x = ufl.SpatialCoordinate(cubeMesh)
ds = ufl.Measure("ds", domain=cubeMesh, subdomain_data=boundaryPanelMeshTagsObj)    

# defining function space for the temperature field (nodes at the corners of the tetrahedra)
nodeSpace = fem.FunctionSpace(cubeMesh, ("Lagrange", 1))

# defining function space for the material properties (nodes at the centre of the tetrahedra)
tetraCentreSpace = fem.FunctionSpace(cubeMesh, ("DG", 0))

# defining function space for the boundary conditions (boundary condition is the heat flux at the centroids of the boundary panels)
# TODO: define correct function space
boundarySpace = fem.FunctionSpace(cubeMesh, ("DG", 0))


# defining the solution at previous time step
u_lastStep = fem.Function(nodeSpace)
u_lastStep.name = "u_lastStep"

# defining the solution at current time step
u_h = fem.Function(nodeSpace)
u_h.name = "u_h"

# apply initial condition
u_h.x.array[:] = np.ones_like(u_h.x.array) * 273
u_lastStep.x.array[:] = np.ones_like(u_h.x.array) * 273

# defining material property field (alpha)
alpha_field = fem.Function(tetraCentreSpace)
alpha_field.name = "alpha_field"
# for this example, we will set the alpha field to a constant value (in reality a function of the material properties and temperature field)
alpha_field.x.array[:] = np.ones_like(alpha_field.x.array) * 0.00018291

# defining material property field (k)
k_field = fem.Function(tetraCentreSpace)
k_field.name = "k_field"
# for this example, we will set the k field to a constant value (in reality a function of the material properties and temperature field)
k_field.x.array[:] = np.ones_like(k_field.x.array) * 413


# defining the source term field (f)
f_field = fem.Function(tetraCentreSpace)
f_field.name = "f_field"
# for this example, we will set the source term to zero
f_field.x.array[:] = np.zeros_like(f_field.x.array)

# defining the boundary condition field (q)
q_field = fem.Function(boundarySpace)
q_field.name = "q_field"
# specifying a random heat flux at the boundary -> in reality this will be an input from a different code module (trajectory propagator)
# TODO: use correct function space for this boundary condition (only the heat flux at the centroids of the boundary panels is known)
q_field.x.array[:] = np.random.uniform(5000000, 8000000, q_field.x.array.shape[0])

# defining the time step
dt = 0.1
t_eval = 0
no_steps = 10

# setup variational problem
u, v = ufl.TrialFunction(nodeSpace), ufl.TestFunction(nodeSpace)

a =  u * v * ufl.dx + dt * alpha_field * ufl.dot( ufl.grad(u), ufl.grad(v)) * ufl.dx
L = (u_lastStep + dt * f_field) * v * ufl.dx

# applying neuamnn boundary condition
L -= dt * alpha_field * -1/k_field * q_field * v * ds


bilinear_form = fem.form(a)
linear_form = fem.form(L)

A = assemble_matrix(bilinear_form)
A.assemble()
        
b = create_vector(linear_form)

solver = PETSc.KSP().create(cubeMesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)



# setup file to write solution to
xdmf = io.XDMFFile(cubeMesh.comm, "diffusion.xdmf", "w")
xdmf.write_mesh(cubeMesh)
xdmf.write_function(u_h, t_eval)    
pyvista.start_xvfb()
grid = pyvista.UnstructuredGrid(*plot.vtk_mesh(nodeSpace))

plotter = pyvista.Plotter()
plotter.open_gif("u_time.gif", fps=10)            
    
grid.point_data["u_h"] = u_h.x.array
warped = grid.warp_by_scalar("u_h", factor=1)        

viridis = mpl.colormaps.get_cmap("jet").resampled(50)
sargs = dict(title_font_size=25, label_font_size=20, fmt="%.2e", color="black",
            position_x=0.1, position_y=0.8, width=0.8, height=0.1)

renderer = plotter.add_mesh(warped, show_edges=True, lighting=False,
                            cmap=viridis, scalar_bar_args=sargs,
                            clim=[0, 2*max(u_h.x.array)])



for i in range(no_steps):
    print(f"Time step {i} / {no_steps} completed \t\t", end="\r")
    
    t_eval += dt
    
    with b.localForm() as loc_b:
        loc_b.set(0)
    assemble_vector(b, linear_form)
    
    
    # Solve linear problem
    solver.solve(b, u_h.vector)
    u_h.x.scatter_forward()

    # Update solution at previous time step (u_n)
    u_lastStep.x.array[:] = u_h.x.array
        
    
    # Write solution to file
    xdmf.write_function(u_h, t_eval)
    # Update plot
    warped = grid.warp_by_scalar("u_h", factor=1)
    plotter.update_scalars(u_h.x.array, render=False)
    plotter.write_frame()
    
plotter.close()
xdmf.close()

I hope this is all clear, if not I’m happy to elaborate more

2

A simple option could be to define a DG 0 field. Note that:

  • it would not be defined only on the boundary, but also in the interior. However, since you never use the value of q in the interior, you could simply leave the values there zero (= the default value)
  • it would not be defined only on the centroid of a cell, but in the entire cell, with constant value in the cell. Still, you’d need to extend the values you get from your engineering model anyways, and extending it in a piecewise constant manner is the simplest possible choice.
1 Like