Non-Convergence in Heat Equation Solver with Temperature-Dependent Cp

General
1

I am solving the heat equation in FEniCS with a temperature-dependent specific heat capacity๐ถ๐‘(๐‘‡). Below is the function I use to define ๐ถ๐‘(๐‘‡) to consider latent heat:

 def Cp(T):
    dependent = True
    if dependent: 
        """Temperature-dependent specific heat capacity (Cp)."""
        # Temperature thresholds
        # Add temperature bound at the start
        T = fe.conditional(fe.lt(T, fe.Constant(298)), fe.Constant(298), T)
        T_solidus = 1400
        T_liquidus =1600
        Cp_solidus = fe.conditional(fe.lt(T, T_liquidus),fe.conditional(fe.lt(T, fe.Constant(298)), .44 ,T*.000173 + .44),.65)
        # Cp_solidus = .5
        Cp_liquid = 0.72
        # Fraction of solid and liquid
        
        f_l = fe.conditional(fe.lt(T, T_solidus), 0,
                            fe.conditional(fe.lt(T, T_liquidus),
                                            (T - T_solidus) / (T_liquidus - T_solidus), 1))
        f_s = 1-f_l
        # Weighted average with latent heat contribution
        Cp_eff = f_s * 0.637994 + f_l *  Cp_liquid
        Cp_eff += (270 / (T_liquidus - T_solidus))  # Add latent heat contribution
        Cp_transition = Cp_eff
        return fe.conditional(
                fe.lt(T, T_solidus),  # Below solidus
                Cp_solidus,
                fe.conditional(
                    fe.lt(T, T_liquidus),  # Between solidus and melting
                    Cp_transition,
                Cp_liquid,  # Above melting point
                ),
            )

The weak form of my equation is as follows:

    theta = 0.5
    u_rhs = theta*u_pre + (1 - theta)*u_crt
    residual = rho*Cp(u_rhs)*(u_crt - u_pre) * v * fe.dx + dt * fe.dot(k* fe.grad(u_rhs), fe.grad(v)) * fe.dx \
                - dt*q_top * v * ds(2) 

        solver_parameters = {'newton_solver': {'maximum_iterations': 20, 'linear_solver': 'mumps'}}
        fe.solve(residual == 0, u_crt, bcs, solver_parameters=solver_parameters)

However, the solver fails to converge. Could this issue be caused by sharp changes in ๐ถ๐‘due to the conditional statements, or is there another potential source of the problem? What would be the best way to address this?

2

Skimming your code it looks like a latent heat problem. If this is the case, these types of coefficients are not straightforward to implement with a โ€œtemperatureโ€ formulation as you have here.

Options which are typically explored:

  1. Smoothing of the coefficients to avoid the degeneracy at the temperature of phase change.
  2. Enthalpy formulation of the underlying model.
  3. A custom solver solving for the latent heat as a problem defined on the quadrature points, e.g., Transient heat equation with phase change. There may be examples from others of which Iโ€™m unaware.