A Coupled PDE-ODE Model for Nonlinear Transient Heat Transfer with Convection Heating at the Boundary: Numerical Solution by Implicit Time Discretization and Sequential Decoupling

: This article considers heat transfer in a solid body with temperature-dependent thermal conductivity that is in contact with a tank ﬁlled with liquid. The liquid in the tank is heated by hot liquid entering the tank through a pipe. Liquid at a lower temperature leaves the tank through another pipe. We propose a one-dimensional mathematical model that consists of a nonlinear PDE for the temperature along the solid body, coupled to a linear ODE for the temperature in the tank, the boundary and the initial conditions. All equations are converted into a dimensionless form reducing the input parameters to three dimensionless numbers and a dimensionless function. A steady-state analysis is performed. To solve the transient problem, a nontrivial numerical approach is proposed whereby the differential equations are ﬁrst discretized in time. This reduces the problem to a sequence of nonlinear two-point boundary value problems (TPBVP) and a sequence of linear algebraic equations coupled to it. We show that knowing the temperature in the system at time level n − 1 allows us to decouple the TPBVP and the corresponding algebraic equation at time level n. Thus, starting from the initial conditions, the equations are decoupled and solved sequentially. The TPBVPs are solved by FDM with the Newtonian method.


Introduction
This work investigates transient heat transfer in a solid body [1] with temperaturedependent thermal conductivity.The body is in thermal contact with a tank filled with a liquid at temperature T. Through one pipe, liquid at a higher temperature is entering the tank.Through another pipe, liquid at temperature T is leaving the tank.Thus, the liquid in the tank is being heated.Systems like this are encountered in engineering and chemical engineering equipment, e.g., in hydronic heating devices and continuous stirredtank reactors (CSTR) [2,3].Due to the symmetry of the considered system, the temperature u in the solid body depends on one spatial variable only; hence, the problem is onedimensional.Since the thermal conductivity of the sloid depends on u, the heat equation for the solid body is a nonlinear partial differential equation (PDE).Due to its importance, the equation has been studied extensively.Various numerical methods for solving the equation have been developed [4][5][6][7].For particular cases, analytical techniques have been proposed [8][9][10].At the liquid-solid contact surface, we apply a convection boundary condition [11].This condition connects the temperature T in the bulk of the liquid with the temperature u and the temperature gradient ∂u/∂x at the surface.In this work, we assume that the liquid in the tank is well stirred so that the temperature in the liquid remains essentially uniform in space throughout the process (perfect mixing, ideal CSTR [2,3]).Consequently, the energy-balance equation of the tank leads to a first-order linear ordinary differential equation (ODE) for the temperature in the liquid T(t), where t is the time.The equation contains the unknown temperature u at the liquid-solid contact surface.Thus, we get a system with a coupled PDE and ODE.Coupled PDE-ODE systems are important in science and engineering and are often studied in connection with control problems of electromagnetic, mechanical, and chemical-reaction coupling [12][13][14][15][16][17].Examples of systems modeled by coupled PDEs and ODEs can be found in [18][19][20].In general, decoupling the PDE and the ODE is not a trivial task, and if possible, leads to a hard-to-handle integrodifferential system [12].In the literature, there are mainly two strategies for numerical solutions of nonlinear PDE-ODE systems.In the first strategy, called coupled, the PDE and the ODE are solved together.Broadly speaking, there are two ways to do this.We can discretize the whole system in the considered space-time domain, e.g., by using the finite difference method (FDM), and then solve the obtained nonlinear algebraic system using an iterative method, e.g., the Newtonian method.Alternatively, we can first linearize the whole system and then solve the obtained sequence of linear PDE-ODE subsystems, e.g., by FDM or the finite element method (FEM).Examples of techniques associated with the first strategy include the nonlinear multigrid method [7,21], the full approximation scheme (FAS) [22], the Newton multigrid method [23,24], and Newton-type linearizationbased methods [25].The second strategy, called decoupled, in the context of the considered problem, is the following.A guess is made for the temperature u on the contact surface as a function of time.Using this guess, we can solve the ODE, obtaining the temperature in the liquid T(t).By substituting T(t) into the boundary condition (BC), we can solve the PDE, subject to the BC, separately from the ODE.The solution provides us with a new function for the temperature u at the surface.The procedure is repeated until convergence.In fact, the method is a fixed-point iteration method.Examples of problems solved by the second strategy can be found in [26][27][28].This method is not recommended for strongly coupled systems.
In this work, we propose a new method to solve the coupled nonlinear PDE-ODE system that does not fall into either of the described strategies.First, analogously to the technique used in [29][30][31], the system is discretized in time using an implicit finite difference scheme.This results in a fully discretized ODE and a semi-discretized PDE.The equations discretized in this way, together with the boundary conditions, can be converted into a sequence of nonlinear two-point boundary value problems (TPBVP) [32][33][34][35][36][37][38][39] for the unknown temperature in the solid and a sequence of linear algebraic equations for the unknown temperature in the liquid.This, as shown in the article, effectively decouples the PDE and the ODE because knowing the temperature in the liquid and along the solid body at time level n − 1 allows us to decouple and solve the TPBVP and the corresponding algebraic equation at time level n.Thus, by starting from the initial conditions and sequentially decoupling and solving the TPBVPs and the algebraic equations, we obtain the solution of the problem.Note that the method is essentially different from the method of lines (MOL) [11,40], which is a traditional numerical approach to solving parabolic (and other) PDEs, because the MOL first discretizes the PDE in space.This yields a system of first-order ODEs with initial conditions, i.e., an initial value problem (IVP).Then, to solve the IVP, one can use Runge-Kutta methods or other appropriate IVP techniques.Contrary to this approach, we first discretize the equations in time, thereby reducing the problem to a sequence of nonlinear TPBVPs.To solve the TPBVPs, we apply the FDM with the Newtonian method [34][35][36], but many other methods can be used.An advantage of the proposed approach is that at each time level, we only have to solve one nonlinear TPBVP.The Jacobian matrix for the Newtonian method is only N × N, where N is the number of the space mesh-points.The method is unconditionally stable, fast, and easy to implement.This article is structured in the following way.Section 2 presents the physical system and the proposed mathematical model.In Section 3, all equations are converted into a dimensionless form, and the mathematical problem is formulated in terms of the dimensionless equations.A steady-state analysis is performed in Section 4. The proposed numerical method is presented in Section 5.In this section, we convert the PDE-ODE system into a sequence of TPBVPs.How to solve the TPBVPs through the FDM is shown in Section 6.In Section 7, several computer experiments are presented, and the results are discussed.

Physical System and Mathematical Model
The physical system considered in the article is the following.A cylindrical solid body is placed along the x-axis between x = 0 and x = L (Figure 1).The cylinder is laterally insulated.The temperature in the solid body is denoted by u(x, t), where t is the time.Due to the symmetry of the problem, the temperature is a function of one spatial variable only.At x = 0, the body is in thermal contact with a tank of volume V filled with liquid.The area of the contact surface is A. The temperature inside the tank is denoted by T(t).The liquid in the tank is well stirred throughout the process, so the temperature is deemed essentially uniform in the bulk.Through one pipe, a hot liquid at temperature T r is entering the tank at a constant volumetric flow rate Q.Through another pipe, a liquid at temperature T(t) is leaving the tank at the same volumetric flow rate.Initially, the temperature in the body and in the tank is T 0 , where T 0 < T r .At x = L, the temperature is kept constant at T 0 , but other boundary conditions at x = L can also be applied.
Newtonian method [34][35][36], but many other methods can be used.An advantage of the proposed approach is that at each time level, we only have to solve one nonlinear TPBVP.The Jacobian matrix for the Newtonian method is only  × , where  is the number of the space mesh-points.The method is unconditionally stable, fast, and easy to implement.This article is structured in the following way.Section 2 presents the physical system and the proposed mathematical model.In Section 3, all equations are converted into a dimensionless form, and the mathematical problem is formulated in terms of the dimensionless equations.A steady-state analysis is performed in Section 4. The proposed numerical method is presented in Section 5.In this section, we convert the PDE-ODE system into a sequence of TPBVPs.How to solve the TPBVPs through the FDM is shown in Section 6.In Section 7, several computer experiments are presented, and the results are discussed.

Physical System and Mathematical Model
The physical system considered in the article is the following.A cylindrical solid body is placed along the -axis between  = 0 and  =  (Figure 1).The cylinder is laterally insulated.The temperature in the solid body is denoted by (, ), where  is the time.Due to the symmetry of the problem, the temperature is a function of one spatial variable only.At  = 0, the body is in thermal contact with a tank of volume  filled with liquid.The area of the contact surface is .The temperature inside the tank is denoted by ().The liquid in the tank is well stirred throughout the process, so the temperature is deemed essentially uniform in the bulk.Through one pipe, a hot liquid at temperature   is entering the tank at a constant volumetric flow rate .Through another pipe, a liquid at temperature () is leaving the tank at the same volumetric flow rate.Initially, the temperature in the body and in the tank is  0 , where  0 <   .At  = , the temperature is kept constant at  0 , but other boundary conditions at  =  can also be applied.The mathematical model describing the heat transfer in the system consists of the following equations.The heat equation for the solid body [1] is where  is the thermal conductivity of the solid body,  is its density, and   is the specific heat capacity at constant pressure.The left-hand side of Equation ( 1) is the rate at which the energy density in the solid body at position  increases.It should be equal to the negative divergence of the heat flux vector at , which in one dimension is just the right-hand side of Equation (1).Equation ( 1) is a PDE for the unknown temperature (, ).In our model,  depends on ; hence, the PDE is nonlinear.The energy balance equation for the tank [2,3] is The mathematical model describing the heat transfer in the system consists of the following equations.The heat equation for the solid body [1] is where κ is the thermal conductivity of the solid body, ρ is its density, and c p is the specific heat capacity at constant pressure.The left-hand side of Equation ( 1) is the rate at which the energy density in the solid body at position x increases.It should be equal to the negative divergence of the heat flux vector at x, which in one dimension is just the right-hand side of Equation (1).Equation ( 1) is a PDE for the unknown temperature u(x, t).In our model, κ depends on u; hence, the PDE is nonlinear.The energy balance equation for the tank [2,3] is Rate of energy change Rate of energy input − AH(T(t) − u(0, t)) Rate of heat loss at the surface where H is the mean convective heat transfer coefficient, which is considered constant.The density of the liquid is ρ l , and its specific heat capacity at constant pressure is c p,l .This is the well-known lumped capacitance model [1] applied to the liquid in the tank [2].It follows naturally from the perfect mixing assumption [3], i.e., when no temperature gradient exists.The left-hand side of Equation (2) is the rate at which the energy in the tank increases.It equals the rate at which energy is entering the tank with the hot incoming liquid, minus the rate at which energy is leaving the tank with the outgoing liquid, minus the rate at which energy is being transferred from the liquid to the solid body through the contact surface.The transport of energy in the pipes is due to bulk motion of the fluid.This transport mechanism is called advection.The transport of energy from the liquid to the solid body is due to convection (convective heat transfer).Equation ( 2) is a linear ODE for the unknown temperature in the tank T(t).This equation is coupled to the PDE in Equation (1) because it contains the temperature at the contact surface u(0, t), which is an output variable (dependent variable) of Equation (1).At x = 0, we apply the following convection boundary condition [1,11]: since κ depends on u, this condition is nonlinear.Note that, since Equation (3) contains the unknown temperature in the tank T(t), we cannot solve the PDE, Equation (1), subject to the boundary condition, Equation (3), separately from the ODE, Equation (2).Hence, the PDE and the ODE are coupled.The initial conditions are at the right boundary, x = L, the temperature is kept constant at T 0 ; i.e., the condition in Equation ( 5) is a Dirichlet boundary condition.Other boundary conditions, e.g., a fixed flux (Neumann-like) or a convection condition (Robin-like), can also be applied at x = L.
In the next section, in order to reduce the number of input parameters, we converted the equations into a dimensionless form.Then, the mathematical problem is formulated in terms of the dimensionless equations.

Dimensionless Equations and Formulation of the Mathematical Problem
The dimensionless temperature in the solid body and in the tank are defined as respectively.Note that u = T 0 and T = T 0 correspond to Θ = 0 and Ψ = 0, respectively; and u = T r and T = T r to Θ = 1 and then, using Equations ( 6)- (8), Equations ( 1)-(3) become −λ(Θ(0, t)) ∂Θ(x, t) Axioms 2023, 12, 323 5 of 22 where both t a and t c have the dimension of time.The rate at which the temperature in the tank Ψ increases equals the difference of the two terms in Equation (10).The first one is due to energy transport by the liquid moving through the pipes (advection) and is inversely proportional to t a .Hence, t a is a characteristic time related to advection.The second term is due to convective heat transfer between the liquid and the solid body through the contact surface at x = 0 and is inversely proportional to t c .Hence, t c is a characteristic time related to convection.The right boundary condition is the requirements that we impose on the thermal conductivity function λ(Θ) is that it is positive and twice continuously differentiable.Let the thermal conductivity at Θ = 0 be λ 0 .
Formally, the solution of the ODE Equation ( 20), subject to the first initial condition in Equation ( 23), can be written as that Equation ( 26) is a solution to Equation ( 20) can easily be checked by direct substation.The problem with Equation ( 26) is that the function θ(0, s) is unknown; hence, the integration cannot be performed.Substituting Equation ( 26) into the boundary condition Equation ( 21) decouples the PDE, Equation ( 19), subject to boundary conditions, Equations ( 21) and ( 22), from the ODE, Equation ( 20), but this leads to an integro-differential system which is hard to handle.In Section 5, we propose a numerical approach for solving the system of Equations ( 19)-( 23).

Steady-State Analysis
Let the steady-state temperature in the tank be ψ s , and in the solid body θ s (ξ): since ψ s is independent of time, Equation ( 20) becomes 1 taking into account that τ a /τ c = t a /t c , the equation can be written as since the function θ s (ξ) is independent of time, the heat equation for the solid body Equation ( 19) becomes let us consider the case g(θ) ≡ 1.Then, taking into account the second boundary condition, θ s (1) = 0, the solution is the temperature gradient is therefore, the convection boundary condition, Equation ( 21), becomes Equation ( 33), which holds for g(θ) ≡ 1, and Equation ( 29), uniquely determine the steady-state temperature in the tank ψ s and on the liquid-solid surface θ s (0).As an example, in Figure 2, we show the steady-state temperature profile for the case where Bi = 3 and t a /t c = 2; hence, by Formulas ( 29) and (33), ψ s = 2/3 and θ s (0) = 1/2.On the right, the temperature intervals [0, θ s (0)], [θ s (0), ψ s ], and [ψ s , 1] are depicted.The ratio of the length of the lower interval to the length of the middle interval equals the Biot number, Bi, Equation (33)-i.e., three.The ratio of the length of the upper interval to the length of the middle interval equals the ratio t a /t c in Equation ( 29), i.e., two.
Let us consider the case () ≡ 1.Then, taking into account the second boundary condi-tion,   (1) = 0, the solution is The temperature gradient is Therefore, the convection boundary condition, Equation ( 21), becomes Equation ( 33), which holds for () ≡ 1, and Equation ( 29), uniquely determine the steady-state temperature in the tank   and on the liquid-solid surface   (0).As an example, in Figure 2, we show the steady-state temperature profile for the case where Bi = 3 and   /  = 2; hence, by Formulas ( 29) and ( 33),   = 2/3 and   (0) = 1/2.On the right, the temperature intervals [0,   (0)], [  (0),   ], and [  , 1] are depicted.The ratio of the length of the lower interval to the length of the middle interval equals the Biot number, Bi, Equation (33)-i.e., three.The ratio of the length of the upper interval to the length of the middle interval equals the ratio   /  in Equation ( 29), i.e., two.

Implicit Time Discretization and Sequential Decoupling
We introduce a time-mesh   = ∆,  = 0,1,2, … , where ∆ is the time-step, and discretize the PDE, Equation (19), and the ODE, Equation ( 20 where   and   are defined in Equation ( 25).The proposed discretization scheme is implicit.This ensures unconditional stability of the numerical method.Let  −1 and the function  −1 () be some suitable approximations for ( −1 ) and (,  −1 ).By dropping the (∆) terms in Equations ( 34) and ( 35) and supplementing the two equations with the boundary conditions Equations ( 21) and ( 22) at  =   , we obtain
In the literature, there are many methods for numerical solutions of nonlinear TPB-VPs [32][33][34][35][36][37][38][39].Some of them linearize the ODE, e.g., by quasi-linearization [37] or Picard linearization [35], and then solve the sequence of linear sub-problems by some linear method, e.g., the collocation method [33] or the linear shooting-method [34].Other methods rely on replacing the TPBVP with a sequence of initial value problems, e.g., the shooting by Newtonian method (nonlinear shooting method) [38] and the shooting-projection method [39].The finite difference method (FDM) [34][35][36] is another important method whereby the ODE is first discretized using finite differences, and then the obtained nonlinear system is solved iteratively, e.g., by the Newtonian method.In this paper, we adopt the FDM because it is reliable, it is easy to implement, and it can handle a large variety of boundary conditions, including the nonlinear condition of Equation (43).Not all available nonlinear methods are appropriate for problems such as Equations ( 42)-( 44), so as explained in [29], caution should be exercised when choosing the particular numerical technique.

Solving the TPBVP via the FDM
We are going to solve the TPBVP in Equations ( 42)-( 44) by the FDM with the Newton iteration method for the solution of the arising nonlinear algebraic systems.For the Newtonian method, we will need the partial derivatives of the function f in Equation ( 45), so we calculate them here.The partial derivative of f with respect to the first variable θ n is the partial derivative of f with respect to the second variable dθ n /dξ is

Discretizing the TPBVP with the FDM
We introduce the space mesh: using the central difference approximations for dθ n (ξ)/dξ and d 2 θ n (ξ)/dξ 2 , the ODE of Equation ( 42) is discretized on the mesh in Equation (48): where in Equation (49) the unknowns θ n,i approximate θ n (ξ i ) with approximation error O (∆ξ) 2 .

Solving the Nonlinear System via the Newton Method
We introduce the vector F n = [F n,1 , F n,2 , . . . ,F n,N ] T with components using Equations ( 54) and ( 55), the nonlinear system in Equations ( 49), (51), and (53) can be written as where n be an approximation of θ n .By expanding where n is the Jacobian of F n with respect to θ n evaluated at θ by replacing in Equation (57) θ n by θ (k+1) n and dropping the higher order terms, we obtain where δθ this is the Newtonian iteration formula.In Equation (59), θ , k = 0, 1, . .., using Equation (59).If the sequence is convergent, then the limiting vector θ n = limθ (k) n , k → ∞ is a solution to the system in Equations ( 49), (51), and (53).In practice, the iteration process is terminated when δθ < , where • is some suitably chosen norm, e.g., the maximum norm, and is some small number.This inequality is called a stopping criterion.
The vector θ (k+1) n is taken as an approximate solution to the system.

Calculating the Elements of the Jacobian
Using Equation (54), we get the nonzero elements of rows i = 2, 3 . . ., N − 1 of the Jacobian equation-Equation (58): in the above formulas, where the first row of the Jacobian equation, Equation (58), concerns the first boundary condition.
Using the first equation in Equation (55), we obtain the nonzero elements of the first row: the last row of the Jacobian concerns the second boundary condition.Using the second equation in Equation (55), we obtain the only nonzero element of the last row: (71)

Computer Experiments and Discussion
In this section, we apply the proposed numerical approach to solve several examples and discuss the obtained results.A MATLAB implementation of the numerical method is given in Appendix A.

Example 1
In this example, we consider dimensionless thermal conductivity of the form problems with exponential dependence of the thermal conductivity on the temperature can be found in [43,44].Such dependence occurs in some materials, e.g., in silicon [43].
In the first experiment, we used Bi = 1 and where τ a and τ c are the dimensionless characteristic times defined in Equation ( 25).The considered time interval is 0, τ f , where τ f = 5.The time interval is discretized by M = 51 mesh-points; hence, the time-step is ∆τ = 0.1.The space-step is ∆ξ = 0.02.At each time level, the Newton iteration is ended when the maximum norm of the difference between two successive approximate solutions has become less than = 10 −6 .Figure 3 shows the transient solution for τ a = 1, τ c = 1, and µ = −2, 0, 2. The results were obtained by the MATLAB program given in Appendix A.
To interpret the results, let us consider the steady-state temperature profiles.At steady state, the heat flux −g(θ s )dθ s /dξ should be the same at any ξ ∈ [0, 1] (see Equation ( 30)).Hence, for a thermal conductivity value that decreases with temperature, we expect the magnitude of the temperature gradient to be high in regions where the temperature is high.This compensates for the low thermal conductivity in these regions.In regions where the temperature is low, we expect low values for the temperature gradient.Conversely, for a thermal conductivity that increases with temperature, the magnitude of the temperature gradient should be high in regions with low temperatures, and vice versa.When the thermal conductivity does not depend on the temperature, we expect a constant temperature gradient.Thus, since for µ = −2, the thermal conductivity equation, Equation (72), is a decreasing function of temperature, we get a convex steady-state temperature profile.For µ = 0, since Equation ( 72) is constant, the steady-state temperature profile is a straight line.For µ = 2, since Equation ( 72) is an increasing function, the steady-state profile is concave.To interpret the results, let us consider the steady-state temperature profiles.At steady state, the heat flux −(  )  / should be the same at any  ∈ [0,1] (see Equation ( 30)).Hence, for a thermal conductivity value that decreases with temperature, we expect the magnitude of the temperature gradient to be high in regions where the temperature is high.This compensates for the low thermal conductivity in these regions.In regions where the temperature is low, we expect low values for the temperature gradient.Conversely, for a thermal conductivity that increases with temperature, the magnitude of the temperature gradient should be high in regions with low temperatures, and vice versa.When the thermal conductivity does not depend on the temperature, we expect a constant temperature gradient.Thus, since for  = −2, the thermal conductivity equation, Equation (72), is a decreasing function of temperature, we get a convex steady-state temperature profile.For  = 0, since Equation ( 72) is constant, the steady-state temperature profile is a straight line.For  = 2, since Equation ( 72) is an increasing function, the steadystate profile is concave.
Next, we investigate how the temperature in the tank and on the liquid-solid contact surface changes with the time  for Bi = 1, a fixed ratio of   /  = 1, and variable values of   .Results are shown in Figure 4.As can be seen in the figure, although the values of the steady-state temperatures for given () are uniquely determined by the Biot number and the ratio of the characteristic times, the individual values of   and   are important and control the rate and manner by which the steady-state is approached.The larger the value of   , the slower the approach.To interpret the results, let us consider the steady-state temperature profiles.At steady state, the heat flux −(  )  / should be the same at any  ∈ [0,1] (see Equation ( 30)).Hence, for a thermal conductivity value that decreases with temperature, we expect the magnitude of the temperature gradient to be high in regions where the temperature is high.This compensates for the low thermal conductivity in these regions.In regions where the temperature is low, we expect low values for the temperature gradient.Conversely, for a thermal conductivity that increases with temperature, the magnitude of the temperature gradient should be high in regions with low temperatures, and vice versa.When the thermal conductivity does not depend on the temperature, we expect a constant temperature gradient.Thus, since for  = −2, the thermal conductivity equation, Equation (72), is a decreasing function of temperature, we get a convex steady-state temperature profile.For  = 0, since Equation ( 72) is constant, the steady-state temperature profile is a straight line.For  = 2, since Equation ( 72) is an increasing function, the steadystate profile is concave.
Next, we investigate how the temperature in the tank and on the liquid-solid contact surface changes with the time  for Bi = 1, a fixed ratio of   /  = 1, and variable values of   .Results are shown in Figure 4.As can be seen in the figure, although the values of the steady-state temperatures for given () are uniquely determined by the Biot number and the ratio of the characteristic times, the individual values of   and   are important and control the rate and manner by which the steady-state is approached.The larger the value of   , the slower the approach.Finally, we present the dimensionless heat flux −g(θ)∂θ/∂ξ through the liquid-solid contact surface, i.e., at ξ = 0, and the dimensionless heat flux through the right boundary, i.e., at ξ = 1, as a function of the time τ for Bi = 0.1, 1, 10, and τ a = τ c = 0.5.To increase the accuracy, we have refined the mesh by using N = 501 space points and M = 501 time points.The results are shown in Figure 5.
Note that, in order to get the heat flux, we have to multiply the dimensionless heat flux by H(T r − T 0 )/Bi.Hence, in Figure 5, it is misleading to compare directly the magnitudes of the fluxes for different Biot numbers.As the results show, the heat flux through the liquidsolid contact surface initially increases quickly, and at about τ = 3 practically "reaches" the steady-state.The corresponding heat flux through the right boundary initially increases very slowly but soon catches up and "reaches", with some delay, the same steady-state value.Although not clearly visible in the figure, close inspection and quantitative results show that the time derivative of the heat flux through the right boundary at τ = 0 is practically zero.For a fixed Biot number, the steady-state value of the heat flux increases with µ.The increase for Bi = 0.1, 1 is almost negligible, but the increase for Bi = 10 is considerable.As can be seen in the figure, for Bi = 10, the steady-state flux for µ = 2 is more than three times greater than the steady-state flux for µ = −2.We can say that the fluxes are more strongly influenced by the nonlinearity parameter µ for high Biot numbers.An interesting phenomenon can be observed for Bi = 10 and µ = −2.Instead of monotonically increasing with time, the heat flux through the liquid-solid contact surface first reaches a maximum and then decreases towards the steady-state value.The phenomenon is not exclusive for high values of Bi and low values of µ, but depending on the values of τ a and τ c , can be observed for other values of Bi and µ.To investigate the phenomenon, we have conducted an experiment for Bi = 1, a fixed ratio of τ a /τ c = 1, and different values of τ a .The results are shown in Figure 6.Finally, we present the dimensionless heat flux −()/ through the liquidsolid contact surface, i.e., at  = 0 , and the dimensionless heat flux through the right boundary, i.e., at  = 1, as a function of the time  for Bi = 0.1, 1, 10, and   =   = 0.5.To increase the accuracy, we have refined the mesh by using  = 501 space points and  = 501 time points.The results are shown in Figure 5.Note that, in order to get the heat flux, we have to multiply the dimensionless heat flux by (  −  0 )/Bi.Hence, in Figure 5, it is misleading to compare directly the magnitudes of the fluxes for different Biot numbers.As the results show, the heat flux through the liquid-solid contact surface initially increases quickly, and at about  = 3 practically "reaches" the steady-state.The corresponding heat flux through the right boundary initially increases very slowly but soon catches up and "reaches", with some delay, the same steady-state value.Although not clearly visible in the figure, close inspection and quantitative results show that the time derivative of the heat flux through the right boundary at  = 0 is practically zero.For a fixed Biot number, the steady-state value of the heat flux increases with  .The increase for Bi = 0.1, 1 is almost negligible, but the increase for Bi = 10 is considerable.As can be seen in the figure, for Bi = 10, the steady-state flux for  = 2 is more than three times greater than the steady-state flux for  = −2.We can say that the fluxes are more strongly influenced by the nonlinearity parameter  for high Biot numbers.An interesting phenomenon can be observed for Bi = 10 and  = −2.Instead of monotonically increasing with time, the heat flux through the liquid-solid contact surface first reaches a maximum and then decreases towards the steady-state value.The phenomenon is not exclusive for high values of Bi and low values of , but depending on the values of   and   , can be observed for other values of Bi and .To investigate the phenomenon, we have conducted an experiment for Bi = 1, a fixed ratio of   /  = 1, and different values of   .The results are shown in Figure 6.Note that, in order to get the heat flux, we have to multiply the dimensionless heat flux by (  −  0 )/Bi.Hence, in Figure 5, it is misleading to compare directly the magnitudes of the fluxes for different Biot numbers.As the results show, the heat flux through the liquid-solid contact surface initially increases quickly, and at about  = 3 practically "reaches" the steady-state.The corresponding heat flux through the right boundary initially increases very slowly but soon catches up and "reaches", with some delay, the same steady-state value.Although not clearly visible in the figure, close inspection and quantitative results show that the time derivative of the heat flux through the right boundary at  = 0 is practically zero.For a fixed Biot number, the steady-state value of the heat flux increases with  .The increase for Bi = 0.1, 1 is almost negligible, but the increase for Bi = 10 is considerable.As can be seen in the figure, for Bi = 10, the steady-state flux for  = 2 is more than three times greater than the steady-state flux for  = −2.We can say that the fluxes are more strongly influenced by the nonlinearity parameter  for high Biot numbers.An interesting phenomenon can be observed for Bi = 10 and  = −2.Instead of monotonically increasing with time, the heat flux through the liquid-solid contact surface first reaches a maximum and then decreases towards the steady-state value.The phenomenon is not exclusive for high values of Bi and low values of , but depending on the values of   and   , can be observed for other values of Bi and .To investigate the phenomenon, we have conducted an experiment for Bi = 1, a fixed ratio of   /  = 1, and different values of   .The results are shown in Figure 6.The results indicate that decreasing the value of τ a while keeping the Biot number and the ratio τ a /τ c fixed, results in a higher and sharper local maximum (peak) in the heat flux through the liquid-solid contact surface.At the same time, the peak shifts towards smaller times.Note that decreasing the characteristic time τ a means that the rate of the energy transfer by advection is increased.This can be achieved by increasing the volumetric flow rate Q.Decreasing τ a while keeping the ratio τ a /τ c fixed means that the characteristic time τ c must also be decreased; i.e., the rate of convective heat transfer must be increased.This can be achieved by increasing H, e.g., by increasing the intensity of the stirring.

Example 2
In this example, we consider the transient behavior of the system for the case the considered time interval is τ ∈ [0, 5], and the time-discretization step ∆τ = 0.1.The final state reached in each experiment is very close to the steady state, so we can investigate the relationship between the dimensionless numbers Bi, t a /t c and the final temperature reached in the tank and on the contact surface, and compare the results with those predicted in Section 4. In the experiments, we used three different Biot numbers: for each Biot number, we have used: (i) τ a = 1/2, τ c = 2 (t a /t c = 1/4), (ii) τ a = 1, τ c = 1 (t a /t c = 1), and (iii) τ a = 2, τ c = 1/2 (t a /t c = 4).Each graph in Figure 7 shows the time evolution of the temperature profile in the solid body (color) and in the tank (black) for a given Biot number Bi and ratio t a /t c .As the results indicate, the temperature in the tank reaches high values when t a /t c is small.This is to be expected since large values of t c , relative to t a , imply a low rate of convective heat transport.Thus, the energy in the tank cannot be transferred fast enough from the tank to the solid body, and the liquid in the tank gets heated considerably.For a fixed Biot number, the lower the value of t a /t c , the higher the steady-state temperature in the tank.As usual, low values of Bi imply a low magnitude of the temperature gradient in the solid body, and high values of Bi, high magnitude of the temperature gradient.Finally, from Figure 7, one can observe that the steady-state temperature in the tank, ψ s , and at the surface, θ s (0), satisfy, as discussed in Section 4, Equations ( 29) and (33).To investigate the convergence of the method, we obtained approximations of the temperature in the tank   and on the liquid-solid surface   | =0 at time  = 1 for where   is the number of space mesh points.The number of time mesh points used in the experiment is  = 10001.The Biot number is Bi = 1, and the characteristic times are are   = 1 and   = 1.Since the exact solution is not available, we assumed that for  = 9, the numerical results are sufficiently close to the exact values, and calculated the errors as To investigate the convergence of the method, we obtained approximations of the temperature in the tank ψ l and on the liquid-solid surface θ l | ξ=0 at time τ = 1 for where N l is the number of space mesh points.The number of time mesh points used in the experiment is M = 10, 001.The Biot number is Bi = 1, and the characteristic times are are τ a = 1 and τ c = 1.Since the exact solution is not available, we assumed that for l = 9, the numerical results are sufficiently close to the exact values, and calculated the errors as The results are shown in Tables 1 and 2. Clearly, since refining the mesh by increasing the number of sub-intervals twice reduces the error four times, the method is second-order accurate with respect to the space step ∆ξ.This means that we can achieve high accuracy using a relatively small number of space mesh-points, and hence, a small Jacobian.When the heat equation is discretized in time using explicit Euler method and in space using central finite differences, we get FDM with the forward time center space (FTCS) scheme [45].For the case considered in this example, Equation (74), the FTCS method is numerically stable if and only if the following condition is satisfied: to compare the stability of our numerical scheme with the FTCS, we solve the problem for Bi Obviously, unlike FTCS, the method is unconditionally stable.When the heat equation is discretized in time using explicit Euler method and in space using central finite differences, we get FDM with the forward time center space (FTCS) scheme [45].For the case considered in this example, Equation (74), the FTCS method is numerically stable if and only if the following condition is satisfied:

Example 3
In this example, we consider three thermal conductivity functions () (Figure 9).The first one is a decreasing function of temperature, the second one is increasing, and the third reaches a maximum at about  = 1/4 and then decreases to a value about 1/5 of its value at  = 0. Accordingly, the obtained steady-state temperature profiles are convex, concave, and a profile with inflexion (Figure 10).In the experiment, we used dimensionless parameters Bi = 4,  = 1/2,  = 2.

Example 3
In this example, we consider three thermal conductivity functions g(θ) (Figure 9).The first one is a decreasing function of temperature, the second one is increasing, and the third reaches a maximum at about θ = 1/4 and then decreases to a value about 1/5 of its value at θ = 0. Accordingly, the obtained steady-state temperature profiles are convex, concave, and a profile with inflexion (Figure 10).In the experiment, we used dimensionless parameters Bi = 4, τ a = 1/2, τ c = 2.

Example 3
In this example, we consider three thermal conductivity functions () (Figure 9).The first one is a decreasing function of temperature, the second one is increasing, and the third reaches a maximum at about  = 1/4 and then decreases to a value about 1/5 of its value at  = 0. Accordingly, the obtained steady-state temperature profiles are convex, concave, and a profile with inflexion (Figure 10).In the experiment, we used dimensionless parameters Bi = 4,   = 1/2,   = 2.

Conclusions
This article considered a PDE-ODE model for nonlinear heat transfer with convection heating at the boundary.For the solution of the coupled PDE-ODE system, a new numerical method was proposed that reduces the system to a sequence of TPBVPs.The method is unconditionally stable, easy to implement, and fast.Computer experiments were performed with various values of the input parameters, demonstrating their influences on the heat transfer process.The steady-state temperatures are uniquely determined by the thermal conductivity function, the Biot number, and the ratio of the two characteristic times.The individual values of the characteristic times, however, control the rate and manner in which the steady state is approached.The smaller the characteristic times, the faster the approach.For small values of the characteristic times, an interesting phenomenon is observed whereby the heat flux through the liquid-solid contact surface reaches a maximum and then decreases towards the steady-state value.Decreasing the characteristic times makes the maximum higher and sharper and shifts it to the smaller times.

Conclusions
This article considered a PDE-ODE model for nonlinear heat transfer with convection heating at the boundary.For the solution of the coupled PDE-ODE system, a new numerical method was proposed that reduces the system to a sequence of TPBVPs.The method is unconditionally stable, easy to implement, and fast.Computer experiments were performed with various values of the input parameters, demonstrating their influences on the heat transfer process.The steady-state temperatures are uniquely determined by the thermal conductivity function, the Biot number, and the ratio of the two characteristic times.The individual values of the characteristic times, however, control the rate and manner in which the steady state is approached.The smaller the characteristic times, the faster the approach.For small values of the characteristic times, an interesting phenomenon is observed whereby the heat flux through the liquid-solid contact surface reaches a maximum and then decreases towards the steady-state value.Decreasing the characteristic times makes the maximum higher and sharper and shifts it to the smaller times.Office-NKFIH, grants K137699 and SNN125119.

Nomenclature
= κ(u), thermal conductivity of the solid body at Θ W/(m•K) λ 0 = λ(0), thermal conductivity of the solid body at Θ = 0 W/(m•K) g(Θ) = λ(Θ)/λ 0 , dimensionless thermal conductivity of the solid body at Θ t a = V/Q, characteristic time related to advection through the pipes s t c = ρ l c p,l V/(H A), characteristic time related to convection in the tank s t = ρc p L 2 /λ 0 , characteristic time related to conduction in the solid body s ξ = x/L, dimensionless position τ = t/t, dimensionless time (Fourier number at θ = 0) θ(ξ, τ) = Θ(x, t), dimensionless temperature in the solid body at ξ and τ ψ(τ) = Ψ(t), dimensionless temperature in the tank at τ g(θ) = g(Θ), dimensionless thermal conductivity of the solid body at θ τ a = t a /t, dimensionless characteristic time related to advection τ c = t c /t, dimensionless characteristic time related to convection Bi = HL/λ 0 , Biot number at reference temperature θ = 0
n is the next, in general better, approximation of θ n .As an initial approximation θ (0) n , we use θ n−1 , i.e., the solution obtained at the previous time level.Starting from θ (0) n , we can find each next approximation θ (k+1) n

Figure 3 .
Figure 3.In the top row, the temperature in the tank ψ(τ) (black) and the temperature in the solid body θ(ξ, τ) (color) for g(θ) = e µθ , µ = −2, 0, 2 are shown.In the bottom row, the corresponding side views are shown, i.e., the temperature profiles at times τ = 0.1, 0.2, . . ., 5.0.Next, we investigate how the temperature in the tank and on the liquid-solid contact surface changes with the time τ for Bi = 1, a fixed ratio of τ a /τ c = 1, and variable values of τ a .Results are shown in Figure 4.As can be seen in the figure, although the values of the steady-state temperatures for given g(θ) are uniquely determined by the Biot number and the ratio of the characteristic times, the individual values of τ a and τ c are important and control the rate and manner by which the steady-state is approached.The larger the value of τ a , the slower the approach.

Figure 5 .
Figure 5.The dimensionless heat flux through the contact surface liquid-solid-body (blue) and through the right boundary (green) as a function of the time  for   =   = 0.5 and three different Biot numbers.The thermal conductivity is () =   ,  = −2, 0, 2.

Figure 6 .Figure 5 .
Figure 6.The dimensionless heat flux through the contact surface liquid-solid-body (blue) and through the right boundary (green) as a function of the time  for Bi = 1, a fixed ratio of   /  =

Figure 5 .
Figure 5.The dimensionless heat flux through the contact surface liquid-solid-body (blue) and through the right boundary (green) as a function of the time  for   =   = 0.5 and three different Biot numbers.The thermal conductivity is () =   ,  = −2, 0, 2.

22 Figure 7 .
Figure 7. Temperature profiles at times  = 0.1, 0.2, … , 5.0 for three different Bio numbers and three different ratios of the characteristic times.The thermal conductivity is () ≡ 1.

22 Figure 10 .
Figure 10.In the top row, the temperature in the tank () (black) and the temperature in the solid body (, ) (color) for various thermal conductivity functions () are shown.In the bottom row, the corresponding temperature profiles at times  = 0.1, 0.2, … , 5.0 are shown.

Figure 10 .
Figure 10.In the top row, the temperature in the tank ψ(τ) (black) and the temperature in the solid body θ(ξ, τ) (color) for various thermal conductivity functions g(θ) are shown.In the bottom row, the corresponding temperature profiles at times τ = 0.1, 0.2, . . ., 5.0 are shown.

Table 1 .
Approximation ψ l of the temperature in the tank at time τ = 1 and the corresponding error E l for different numbers of space mesh points N l .

Table 2 .
Approximation θ l | ξ=0 of the temperature on the liquid-solid surface at time τ = 1 and the corresponding error e l for different number of space mesh-points N l .

Table 2 .
Approximation   | =0 of the temperature on the liquid-solid surface at time  = 1 and the corresponding error   for different number of space mesh-points   .

Table A1 .
Notation and corresponding MATLAB variables.