Implicit time discretization and quasilinearization for numerical solution of transient nonlinear heat transfer problems Open Access

The one-dimensional heat equation (Cannon, 1984; Carslow et al., 1986; Bergman, 2011) with temperature-dependent thermal conductivity is a nonlinear partial differential equation (PDE) (Polyanin and Zaitsev, 2012). The exact solution to this equation has been found only for a limited number of specific cases (Polyanin et al., 2000; Sadighi and Ganji, 2007; Wazwaz, 2001, 2003; Kosov and Semenov, 2019). For this reason, some approximate analytical techniques (Mustafa et al., 2014; Hristov, 2015, 2016; Yu et al., 2024, and the references therein) and a variety of numerical methods (Bartish and Shakhno, 1993; Gürarslan, 2010; Filipov and Faragó, 2018; Zheng et al., 2020; Lienemann et al., 2005; Gürarslan and Sari, 2011; Zhuang et al., 2021; Zen et al., 2024; Feng et al., 2020) have been developed. A typical numerical approach is to fully discretize the heat equation, or a linearized version of it, using the finite difference method (FDM) (Smith, 1985; LeVeque, 2007) or the finite element method (Taler and Ocłoń, 2014), whereby a system of algebraic equations is obtained. Then, numerical methods for algebraic systems such as Gauss–Seidel (Burden et al., 2016) or Newton’s method (Quarteroni et al., 2010) can be used. Efficient solution of systems resulting from discretizing PDEs is provided by the multigrid methods (Brabazon et al., 2014; Brandt, 1977; Briggs et al., 2000; Trottenberg et al., 2001) such as the full approximation scheme (Brabazon et al., 2014; Brandt, 1977) or the Newton-Multigrid (Brabazon et al., 2014; Briggs et al., 2000; Trottenberg et al., 2001). Note that, because the original problem is nonlinear, an iterative method has to be involved at some stage of the numerical solution. Along with other possibilities, Newton and Picard iterations are often used (Zen et al., 2022). An alternative to full discretization of the PDE is a semi-discretization approach, whereby the PDE is discretized only with respect to one independent variable. This converts the PDE into a system, or sequence, of ordinary differential equations (ODEs). Then, one can use numerical methods for ODEs (Ascher and Petzold, 1998).

An important numerical method that relies on semi-discretization is the method of lines (Liskovets, 1965; Zafarullah, 1970; Schiesser, 2012; Marucho and Campo, 2016). It was originally introduced by Liskovets (1965) for PDEs of elliptic, parabolic, and hyperbolic type. Since then, it has been applied in many different contexts and is still actively used today (Marucho and Campo, 2016). The method discretizes the equation in space, whereby a system of first-order ODEs is obtained. The system, together with the initial condition, constitutes an initial value problem (IVP). It can be solved using various numerical approaches for IVP (Ascher and Petzold, 1998), which, of course, may involve further discretization.

Recently, another approach based on semi-discretization has been proposed (Filipov and Faragó, 2018). Initially, using an implicit scheme, the method discretizes the nonlinear heat equation in time only, whereby, adding the boundary conditions, a sequence of nonlinear two-point boundary value problems (TPBVPs) is obtained. Starting from the initial condition, the TPBVPs can be solved successively using a numerical method for nonlinear boundary value problems for ODEs (e.g. Cuomo and Marasco, 2008, and the references therein). In Filipov and Faragó (2018), to solve the nonlinear TPBVPs, the FDM with Newton’s method has been applied. The original work considers only Dirichlet boundary conditions, but the method has also been applied successfully to relaxing boundary conditions (Filipov et al., 2022) and PDE–ODE systems (Filipov et al., 2023).

The work presented in this article generalizes the original method (Filipov and Faragó, 2018) to nonlinear boundary conditions that can depend on time explicitly. Moreover, contrary to the original approach, which applies the FDM directly to the nonlinear TPBVP and then uses the Newton method, we first apply the quasilinearization method (Bellman and Kalaba, 1965; Lakshmikantham et al., 1995; Mohapatra et al., 1997; Ahmad et al., 2001; Khan, 2006) and just then the FDM. As proven in Filipov et al. (2019), these two approaches formally lead to the same result; however, the latter approach has an advantage. Because the quasilinearization method replaces the nonlinear TPBVP by a sequence of linear TPBVPs, one can take advantage of the many available methods for numerical solution of linear problems such as the FDM, the collocation method, and others (Ascher et al., 1995). The proposed method is computationally very efficient. The time complexity of the algorithm is O(MN), where M is the number of time levels and N is the number of space mesh points. Moreover, the method is unconditionally stable, i.e. it remains numerically stable regardless of the chosen values for M and N.

Let u(x, t) be the temperature at position x and time t, and α be the thermal diffusivity of the media:

where ρ is the density, C_p is the heat capacity at constant pressure, and κ is the thermal conductivity. The density and the heat capacity are considered constant, but the thermal conductivity is assumed to depend on the temperature (e.g Yu et al., 2024; Ivanova et al., 2020); hence, α = α(u). Then, the heat equation is:

We assume that α is positive and two times continuously differentiable. Obviously, unless α is constant, the PDE (1) is nonlinear. The equation is considered on the domain D ={(x, t)|a < x < b, t > 0}, i.e. in the space region between x = a and x = b, and for t > 0 (Figure 1).

At the boundaries, we impose boundary conditions of the form:

where g_a and g_b are given – not necessarily linear – functions that have continuous partial derivatives with respect to the first and second argument. The initial condition is:

where u₀ (x) is the given initial temperature profile. Performing the differentiation on the right-hand side of equation (1), we obtain:

where D_uα = dα/du. We introduce the time mesh t_n = nτ, n = 1,2,… Let us consider equation (5) at the time level t = t_n. Using the Taylor expansion u(x, t_n−1) = u(x, t_n) − τ∂_t u(x, t_n) + O(τ²), we replace the time derivative in (5) by the backward finite difference approximation plus O(τ) term:

Therefore, dropping the O(τ) term, we get the following semi-discretization in time of the PDE (1):

where ' denotes differentiation with respect to x. If u_n−1 (x) is a given (known) function, then equation (6) is a second-order nonlinear ODE for the function u_n (x). The time-discretization scheme used in equation (6) leads to the backward Euler method. We choose this method because it is the simplest implicit method, and implicit methods, e.g. backward Euler (Filipov and Faragó, 2018) and Crank–Nicolson (Zen et al., 2024), typically lead to unconditional stability of the whole numerical procedure. In a more compact form, equation (6) can be written as:

where:

Equation (7) and the boundary conditions:

constitute a nonlinear TPBVP for the unknown function u_n (x). In fact, because n = 1,2,…, we have a sequence of nonlinear TPBVPs. Starting from the initial condition u₀ (x), we can solve the TPBVP (7)–(10) for n = 1, i.e. at time level t₁, and get u₁ (x), which is an approximation of u(x, t₁). Then, using u₁ (x), we can solve the TPBVP for n = 2 and get u₂ (x), which is an approximation of u(x, t₂), and so on (Figure 2). Thus, successively, we obtain an approximation u_n (x) of the exact solution at each time level t_n.

Finally, we note that the time-discretization (6) is consistent with the PDE (5) to first order with respect to τ. Hence, the sequence of nonlinear TPBVPs (7)–(10) is consistent with (5) to first order, too. Because we use an implicit method which is typically stable, our approach results in a convergent algorithm.

In this section, we construct an iterative procedure for solving the nonlinear TPBVP (7)–(10) at each time level t_n, n = 1,2,…. The procedure is based on the quasilinearization method. Let $un(k)(x)$ be some given (known) function that approximates the exact solution u_n (x) of the problem (7)–(10). It is called the k-th approximation of the solution. Using $un(k)(x)$⁠, we can expand the right-hand side of equation (7) in Taylor series around the point $(un(k)(x),(un(k)(x))′)$ getting:

where q and p are the partial derivatives of f (8) with respect to u_n and $u′n$⁠, respectively:

Dropping in equation (11), the terms higher than linear and replacing the exact solution of (7)–(10)u_n (x) by its approximation $un(k+1)(x)$⁠, we get the quasilinearization equation:

Equation (14) is a linear ODE for the unknown function $un(k+1)(x)$⁠. Together with the boundary conditions:

it constitutes a linear TPBVP for $un(k+1)(x)$⁠. If the temperature profile at time level n − 1, i.e. u_{n− 1} (x), and $un(k)(x)$ are given (known), we can solve the linear TPBVP (14)–(16) and get the next (k + 1)-th approximation $un(k+1)(x)$⁠. In general, $un(k+1)(x)$ is a better approximation of u_n (x) than $un(k)(x)$⁠. Starting from some zeroth approximation $un(0)(x)$⁠, we can solve the TPBVP (14)–(16) for k = 1,2,… and get a sequence of functions $un(0)(x),un(1)(x),un(2)(x),…$ If the sequence is convergent, then the limiting function $un(x)=limk→∞un(k)(x)$ is the sought solution of (7)–(10). At each time level n as a zeroth approximation $un(0)(x)$⁠, we can take the numerical solution of the nonlinear problem at the previous time level, i.e. u_{n− 1} (x). For time level n = 1, it is just the initial condition u₀ (x). In Figure 3, we have shown an example with time discretization step τ = 1. The conductive media is in the region 0 < x < 1 and is initially at temperature u = 0. At time t = 0, the temperature at x = 0 is suddenly raised to u = 1 and kept fixed at this value for t > 0.

The quasilinearization method is, in fact, a Newton–Raphson method on operator level (Filipov et al., 2019). Its convergence is quadratic. When the zeroth approximation is “sufficiently good,” as in our case, the method converges rapidly (Cuomo and Marasco, 2008). Note that the smaller the time step τ gets, the closer the zeroth approximation gets to the solution (see Figure 3) and the more efficient the quasilinearization method becomes. For the examples considered in Section 7, when the time step is τ = 0.01, it takes only about two to three iterations at each time level to reach the solution with precision 10⁻⁶.

To solve the linear TPBVP (14)–(16), we are going to use the FDM. We introduce the space-mesh:

Let $un,i(k)$ and $un,i(k+1)$ denote approximations of $un(k)(xi)$ and $un(k+1)(xi)$⁠, respectively. For the inner mesh points, the first derivatives of $un(k)(x)$ and $un(k+1)(x)$ at x = x_i are approximated by the central difference approximation:

The error that we make by replacing the derivatives with their approximations (17) is O(h²), i.e. second order in h.

At the left boundary x = a, the first derivatives of $un(k)(x)$ and $un(k+1)(x)$ are approximated by the forward difference approximation with second-order error:

At the right boundary x = b, the first derivatives of $un(k)(x)$ and $un(k+1)(x)$ are approximated by the backward difference approximation with second-order error:

Taking into account (12), (13), and (8), respectively, we introduce the notation:

Using the central difference approximation for the second derivative of $un(k+1)(x)$ at x = x_i, i = 2, 3,…, n − 1, equation (14) is discretized on the mesh, yielding:

The boundary conditions are:

Equations (20)–(22) constitute a system of algebraic equations for the unknown $un,i(k+1),i=1,2,…,N$⁠. They replace the TPBVP (14)–(16). To solve (20), subject to (21)–(22) in linearized form, we introduce the vectors:

where:

Then, according to Theorem 1 in Filipov et al. (2019), the solution of the system is:

where:

and:

is the Jacobian of $Fn(k)$ with respect to $un(k)$⁠. The solution $un(k+1)$ is the next, in general better, approximation of u_n. Starting from the 0-th approximation $un(0)$⁠, i.e. the solution u_n−1 obtained at the previous time level, we can find each next approximation $un(k+1)$⁠, k = 0,1,… using (26). If the sequence is convergent, then the limiting vector $un=limk→∞(un(k))$ is the approximate solution to the nonlinear TPBVP (7), (9)–(10). In practice, the iteration process is terminated when $‖δun(k)‖<ϵ$⁠, where ‖·‖ is some suitably chosen norm, e.g. the maximum norm, and ϵ is some small number. This inequality is called stopping criteria. The vector $un(k+1)$ is taken as approximate solution to the nonlinear TPBVP. As shown in the next section, the Jacobian matrix (28) is tridiagonal or close to tridiagonal for the most common types of boundary conditions. Hence, using the Thomas method and, if necessary, the Sherman–Morrison algorithm, (27) can be computed in only O(N) operations, where N is the number of space mesh points (Filipov et al., 2019). Hence, the time complexity of the whole method is O(MN), where M is the number of time mesh points.

To be able to use iteration (26), we need to have the elements of the Jacobian matrix (28). In this section, we calculate explicitly these elements. The Jacobian matrix is an n × n matrix, where n is the number of space mesh-points. Rows 2 through n – 1 relate to the discretized differential equation, while rows 1 and n relate to the left and right boundary conditions, respectively.

According to equations (23) and (17), the nonzero elements of rows i = 2, 3,…, n − 1 of the Jacobian are:

The first row of the Jacobian concerns the first boundary condition. According to (24) and (18), the nonzero elements of the first row are:

The last row of the Jacobian concerns the second boundary condition. According to equations (25) and (19), the nonzero elements of the last row are:

This section considers several types of boundary conditions that are likely to occur in practice. First, we consider the case where the temperature u at the boundary is specified. This means that the temperature at the boundary is a given function of time. When this function is just a constant, we get the usual Dirichlet boundary condition. If the temperature at the boundary increases or decreases in a monotonic way approaching a certain value as time goes to infinity, we have a relaxing boundary condition (Filipov et al., 2022; Hristov, 2022). If the temperature at the boundary oscillates without approaching a limit, we get an oscillating boundary condition. Examples for these three cases are presented in Section 7.

Let the temperature at the first boundary x = a and the second boundary x = b be u_a (t) and u_b (t), respectively, i.e.:

Moving all terms in (29) and (30) to the left-hand side of the equations and comparing with the boundary conditions (2) and (3), we see that:

Obviously, in the case of specified temperature, the functions g_a and g_b do not depend on the temperate gradient ∂u/∂x and are linear with respect to the temperature u. Hence, in this case, we deal with linear boundary conditions. The derivatives of g_a and g_b with respect to u are: