Using FiPy for internal heat flux boundary conditions

Question

I need some help with setting boundary conditions for a material simulation, which is thermally degrading. The following aspects are of main concern to me at this stage:

• Mesh motion or deformation is substituted by some kind of loose level set method. Therefore a cell can "switch" the domain (from solid to gas phase e.g.), since the material is only degrading from the outside to the inside. This is marked by a level variable. Here, I use internal boundary conditions as described here to describe the moving boundary. This is implemented in the example below.

• I want to set a specified heat flux as boundary condition for the energy equation in form of a temperature equation (see example below). My first attempt would be to use the fixed flux condition described here, but I'm unsure whether this is physically correct. Is it possible to set a specified heat flux boundary also at internal faces, similar to the fixed value/gradient in the link above? How do I describe this in the energy equation? Moreover, there are source terms which do not use existing fipy classes (see example below). Must they be adusted in some way to allow for the internal flux boundary? I was thinking about adding a generic unity coefficient to be set to 0 on the respective faces as shown with the fixed flux boundary in the usage doc.

The minimal working example can be found below. Since it shall only give the frame for my technical questions, it is not physically correct. I'm aware that the units of the energy equation do not match in this example.

from fipy import Grid2D, CellVariable, TransientTerm, DiffusionTerm, Viewer, ImplicitSourceTerm, FaceVariable

# Consider a 2D mesh with initially constant temperature
mesh = Grid2D(nx=50, ny=50)

# Variables are defined below. Main variables are the temperature and the level variable
# Further variables are defined to include a generic term in the equation
temperature = CellVariable(mesh=mesh, name="temperature", value=300.)
pressure = CellVariable(mesh=mesh, name="pressure", value=1.2 * temperature * 330.)
transCoeff = CellVariable(mesh=mesh, name="transient coefficient", value=.01)
diffCoeff = CellVariable(mesh=mesh, name="diffusion coefficient", value=1.)

# Define first mesh row as boundary mask to simulate initial boundary condition
level_variable = CellVariable(mesh=mesh, name="level variable", value=0)
level_variable.setValue(1, where=((mesh.cellCenters[1] > 49) | (mesh.cellCenters[0] < 1)))

# Simple energy equation including a generic convection term and internal boundary terms
largeValue = 1e10
value = 1000.
transCoeff.constrain(0., mesh.exteriorFaces)
diffCoeff.constrain(0., mesh.exteriorFaces)

eqn = (
        TransientTerm(coeff=transCoeff, var=temperature) == DiffusionTerm(coeff=diffCoeff, var=temperature)
        # generic gas convection term
        - (1e2 * pressure.faceGrad).divergence

        # Internal fixed value boundary condition
        # This is to be changed for a "internal fixed flux boundary condition"
        - ImplicitSourceTerm(level_variable * largeValue, var=temperature)
        + level_variable * largeValue * value

        # Fixed Flux boundary condition on exterior faces:
        # + (mesh.facesLeft * exteriorFlux).divergence
)

# Define viewer and run simulation
viewer = Viewer(vars=[temperature, level_variable], cmap="jet")
steps = 200
dt = 1e-2

for i in range(steps):
    # Here, 950K is the criterion for the material to completely decompose
    # After reaching 950K, the level variable is set to 1 and the cell is immediately considered as boundary for 1000 K
    level_variable.setValue(1, where=(temperature >= 950.))
    eqn.solve(var=temperature, dt=dt)
    viewer.plot()

Besides the questions above I would like to understand the physical/mathematical meaning of the internal or fixed flux conditions added within the PDE. Are there any references that can be recommended in this context? I wasn't able to find further infos in the FiPy documentation. Also, general notes on the procedure are very welcome. Physical integrity is as already mentioned not given in the example. Of course, an additional mass conservation (therefore sweeping) and a detailed material model need to be implemented.

Thanks in advance for any help!

Mathematical expression

Consider following PDE. Since the coefficients of the transient terms are outside of the temporal derivative and also CellVariables, the transient term is modelled as follows:

# Transient term:

TransientTerm(var=temperature, coeff=(rho * Cv))
- temperature * (rho * Cv - rho .old * Cv.old)

The mathematical expression herefore is given here. The heat conduction is a straight-forward diffusion term with conductivity as coefficient:

DiffusionTerm(var=temperature, coeff=k)

The last rhs term describes energy change due to some generic phase change. Since there is a temporal change of a CellVariable with is not the var of the equation, I modelled it as follows:

(rho_h * enthalpy_h - rho_h.old * enthalpy_h.old) / dt