Heat transfer equation with temperature-dependent thermal conductivity

The implemented heat transfer equation in FEATool Multiphysics is shown below.

rho_ht*cp_ht*T' - k_ht*(Tx_x + Ty_y) + rho_ht*cp_ht*(u_ht*Tx_t + v_ht*Ty_t) = q_ht

However, this equation is no longer applicable when considering the thermal conductivity as a function of temperature. The equation above becomes:

rho_ht*cp_ht*T' - k_ht*(Tx_x + Ty_y) - (dk_ht/dT)*((Tx_t^2) + (Ty_t^2))  + rho_ht*cp_ht*(u_ht*Tx_t + v_ht*Ty_t) = q_ht.

How can the term dk_ht/dT be included in the heat transfer equation?

The implemented heat transfer formulation should account for non-linear and temperature dependent thermal conductivity. For better or worse the FEATool equation syntax is designed to look like the strong equation while actually being the weak PDE formulation, for example the term "Tx_x" (in contrast to "Txx" without an underscore) is equivalent to "dT/dx * dv/dx" where "v" is a FEM test function (see the table in the section for Custom Equations of the FEATool documentation).

In particular, in the traditional Galerkin FEM derivation the term

\[ \nabla\cdot(k \nabla T) = 0\]

is multiplied by a test function v and integrated

\[ \int_\Omega\nabla\cdot(k \nabla T)\cdot v\ d\Omega = 0\]

after applying partial integration we get

\[ -\int_\Omega k \nabla T\cdot \nabla v\ d\Omega + \int_{\partial\Omega}k\ \hat{n}\cdot \nabla T\ v\ dS = 0 \]

where now the divergence operator has been transferred to the test function space, and the left term above is equivalent to the FEATool syntax, and the right is the ensuing natural Neumann bounary condition.