where \(\rho=\rho(x)\) is mass density, \(c=c(x,T)\) is the specific heat, \(T=T(x,t)\) is the temperature, \(K=k(x,T)\) is the thermal conductivity, \(Q=Q(x,T,t)\) is the calorific value. \(x\) is the position in the modeling domain, \(T\) is the temperature and \(t\) is the time.
The modeling domain is represented by S, and the boundary is represented by \(\varGamma\). When assuming the boundary conditions of either the Dirichlet or Neumann type, those boundary conditions can be mathematically expressed as
where the term \(T_{1}\), \(q\) is already known. \(q\) is the heat flux outflow from the boundary. Three types of heat flux can be considered in WELSIM thermal module.
\[ q=-q_{s}+q_{c}+q_{r} \]
\[ q_{s}=q_{s}(x,t) \]
\[ q_{c}=hc(T-T_{c}) \]
\[ q_{r}=hc(T^{4}-T_{r}^{4}) \]
where \(q_{s}\) is the distributed heat flux, \(q_{c}\) is the heat flux by the convective heat transfer, and \(q_{r}\) is the heat flux by the radiant heat transfer. The other quantities are
\(T_{c}=T_{c}(x,t)\) Convective heat transfer coefficient ambient temperature
\(h_{c}=h_{c}(x,t)\) Convective heat transfer factor
\(T_{r}=T_{r}(x,t)\) Radiant heat transfer coefficient ambient temperature
\(h_{r}=\epsilon\sigma F=h_{r}(x,t)\) Radiant heat transfer factor. \(\epsilon\) is the radiant rate, \(\sigma\) is the Stefan-Boltzmann constant, \(F\) is the shape factor.
Equation \(\eqref{eq:ch5_thermal_gov2}\) is nonlinear and unsteady. When the time is discretized by the backward Euler's rule and the temperature at time t=t_{0} is known, the temperature at t=t_{0+\triangle t} is calculated using the following equation.
Substituting equations \(\eqref{eq:ch5_thermal_gov_disc2}\), \(\eqref{eq:ch5_thermal_gov_disc3}\), and \(\eqref{eq:ch5_thermal_gov_disc4}\) into equation \(\eqref{eq:ch5_thermal_gov_disc1}\) and skipping the high order polynomial terms, we have
Since the implicit time solver is applied in the program, the selection of incremental time \(\triangle t\) is relatively flexible. However, if the magnitude of \(\triangle t\) is too large, the convergence frequency will be decreased in the iterative computation. The program contains automatic incremental functions to monitor the size of the residual vectors during the iterations. As the convergence rate becomes slow, the incremental time \(\triangle t\) is automatically reduced. When the convergence rate becomes high, the program increases the incremental time \(\triangle t\). Doing this automatic scheme can improve the numerical performance and saving computational time.