Numerical Solution of the Heat Transfer Equation Coupled with the Darcy Flow Using the Finite Element Method

The ﬁ nite element approach was utilized in this study to solve numerically the two-dimensional time-dependent heat transfer equation coupled with the Darcy ﬂ ow. The Picard-Lindelöf Theorem was used to prove the existence and uniqueness of the solution. The prior and posterior error estimates are then derived for the numerical scheme. Numerical examples were provided to show the e ﬀ ectiveness of the theoretical results. The essential code development in this study was done using MATLAB computer simulation.


Introduction
Let Ω ⊂ ℝ 2 be an open domain with a bounded boundary that is simply connected in ℝ 2 , with a continuous Lipschitz boundary ∂Ω. The equation of heat transfer coupled with the Darcy flow is based on heat transfer in porous medium [1]. The governing equation for the transfer of heat together with the Darcy flow is stated in a system of partial differential equations as where h and T are the unknowns and q = −K · ∇h is the Darcy velocity. Such a problem is often investigated using the partial differential equation scheme, and exact methods such as variable separation, variable change, and transformations are commonly employed in order to find exact solutions in a particular domain. However, these approaches cannot be used when dealing with nonhomogeneous initial or boundary conditions, nonlinear partial differential equations, and irregular domains or when resources are limited.
As a result, many numerical approaches are utilized to get approximate solutions of equations whose relevance assists in the issues they describe or depict. Many papers have been written on the subject of heat convection in liquid media, in which its movement is described by the Navier-Stokes equations in combination with the heat equation [2][3][4][5][6][7][8]. An alternative coupling Darcy's model with the heat equation with constant viscosity but varying external force as a function of temperature has been studied by [9,10] and discretized with a spectral method. Ganapathy and Mohan [11] studied the mechanics of the steady-state evolution of motion of the fluid and transfer of heat produced by a difference in temperature in geometry under the assumption that the Darcy flow model was valid. They discovered that the motion of the fluid on the free surface behaves similarly to axial flow. He et al. [12] properly built, analyzed, and compared the control equations of porous media under the effects of external energy sources. Analog simulation is done on the model created with COMSOL multiphysical field processing software, taking into consideration the effect of various physical and geometrical characteristics of porous medium on the model. The simulation shows that when external energy sources are present, the porous media are pushed ahead along the transmission direction of external energy in the phase change range on the macro level. Teskeredžić et al. [13] implemented a numerical technique for flow of fluid, transfer of heat, and stress analysis in phase change issues. The approach relies on numerical meshes made up of cells of any polyhedral shape solving integral form equations regulating momentum, mass, and energy balance. Ahmad et al. [14] looked at numerically solving nonlinear differential equations for heat transmission in micropolar fluids across a stretching domain. With proper consideration of micropolar fluid theory, this study delivers realistic and distinct results.
El-Hakiem et al. [15] solved the problem of hydromagnetic dispersion in non-Darcy mixed convection heat and mass transfer from a vertical surface imbedded in porous medium using a similarity solution. Furthermore, Mansour and El-Shaer [16] investigated the effect of a magnetic field on non-Darcy axisymmetric free convection in a power law fluid saturated porous media with variable permeability.
Motivated by the preceding findings, the finite element method is introduced to solve numerically the governing equations (1)- (2). The finite element technique is a computational method for obtaining approximate solutions to partial differential equations (PDE).

Mathematical Model Analysis
2.1. Initial Boundary Condition and Basic Principles. Equations (1)-(2) should be solved by specifying an initial condition: as well as boundary conditions, which could also take the following form: where T 0 , h 0 , p, q, r, and s are assigned functions and f∂Ω D , ∂Ω N g specifies a boundary partition, that is, ∂Ω D ∪ ∂Ω N = ∂ Ω and ∂Ω D ∩ ∂Ω N = ∅. ∂Ω N is the Neumann boundary, and ∂Ω D is the Dirichlet boundary. An integral domain with the solution belonging to the Sobolev space H 1 ðΩÞ and the underlying space L 2 ðΩÞ is utilized, since the problem cannot be expressed in terms of second derivatives. The L 2 ðΩÞ is a space that has the following definition [17].
and the Sobolev space H 1 ðΩÞ is defined as The space for vanishing boundary values is defined as Sobolev's imbedding is used; for all real numbers p ≥ 1, there exist constants S 0 p and S p such that Equation (10) simplifies to Poincaré's inequality when p = 2.

Weak Formulation.
In order to find the numerical solution of the governing equations (1)-(2) using the finite element method, variational formulation is introduced. This is formally proceeding, by multiplying the governing equations (1)-(2) for each t > 0 with a basis function v = vðx, yÞ and integrating on Ω. V = H 1 ðΩÞ is set, and hðtÞ, TðtÞ ∈ V is sought for each t > 0 such that ð : ð :  By putting the Darcy law, that is, q = −K∇h, into (11) and (12), we obtain ð : ð : Then, by using Green's theorem stated in [18], equations (13) and (14) can be ð : ð : In (16), since we considered basis functions to be a linear basis function (e.g., Assume Bringing (19) and (20) into (15) and (16), respectively, we obtain ð : ð : Now considering the Galerkin approximation of the variational formulation (23) and (24) for all t > 0, let V l be the space of vanishing Lagrange finite elements of degree r on the boundary with respect to a mesh of mesh size l and define h l and T l as an approximate solution, h l ðtÞ, T l ðtÞ ∈ V l such that with T l ð0Þ = T 0l and h l ð0Þ = h 0l , where V l ⊂ V is an appropriate finite-dimensional space and T 0l and h 0l are convenient approximations of T 0 and h 0 in the space V l , respectively. The weak formulation (25) and (26) is called semidiscretization of (13) and (14).   Abstract and Applied Analysis porous media using the finite element technique. That is the solution of (26) using the finite element technique. Since H 1 ðΩÞ is detachable, as a result, it has a quantifiable basis fϕ i g i≥1 : Let V l denote the region covered by the first l basis functions, fϕ i g 1≤i≤l : The reduced problem (26) is discretized in V m by the square system of equations. That is, find T l = ∑ 1≤i≤l T i ϕ i ∈ V l , the solution of The Picard-Lindelöf Theorem [19] guarantees the existence and uniqueness of the solution of the variational formulation (27) in which the theorem is stated as follows.
Then, there exists a unique solution T ∈ L 2 ðð0, T f Þ, VÞ of the initial value problem: The detail of the proof is found in [19]. Condition (28) is a term that is frequently used to describe positivity or coercivity, while (29) is called boundedness or continuity.
2.4. Coercivity/Positivity. We need to constrain the flow velocities to guarantee that bðv, vÞ is positive. That is what we require: and assume the bounded velocities kqk ∞ : Using the Poincaré inequality, kvk H 1 0 ≤ C Ω kvk H 1 , where C Ω is some constant [19].

Continuity/Boundedness. The second condition the boundedness follows by applying the Cauchy-Schwartz inequality:
where μ 2 = ðð1/PeÞ + C Ω kqk ∞ Þ > 0: Hence, the continuity condition is satisfied. Therefore, the solution of the weak formulation of our problem exists and is unique.

Approximate Solution.
Let e denote the element number in a region Ω. To provide an algebraic interpretation of (25) and (26), a basis ϕ j is introduced for V l , and it is observed that it suffices that (25) where (i) M e is the 3 × 3 matrix whose elements are given by m e ji = Ð Ω e ϕ j ϕ i dΩ e (ii) A e is the 3 × 3 matrix whose elements are obtained from (19) as a e ji = Ð Ω e K∇ϕ j ∇ϕ i dΩ e (iii) B e is the 3 × 3 matrix whose elements are obtained from (20) as b e ji = Ð Ω e q∇ϕ j ϕ i dΩ e + ð1/PeÞ The matrix A e is called the local stiffness matrix, and M e is called the local mass matrix. Now, by the assembling process, (37) and (38) will be used to get the overall approximate solution and after some algebraic manipulation, which will follow shortly, allow to be put into the following form: in which h = ½h 1 , h 2 , ⋯, h N T , T = ½T 1 , T 2 , ⋯, T N T and the N × N matrix M is the global mass matrix, the N × N matrix A is the global stiffness matrix, and the vectors g and r are the globalized force vectors. Here, N denotes the total number of nodes in the problem as a whole, and the components of h and T are now labeled by their global node numbers.
For the numerical solution of the ODE system (39) and (40), many methods are available from that the backward Euler method is used and is given in the following form: which is first order accurate with respect to Δt = t k+1 − t k . media. Before analyzing a weak formulation (26), its convergence is analyzed, since it is less involved. The key to the analysis of the weak formulation is to compare T l not directly to T, but rather to an appropriate representative w l ∈ C 1 ð½0, T f , V l Þ. For w l , we choose the elliptic projection of T, defined by

Error
From [20], by the finite element method for elliptic problems, we have the L 2 estimate: for some constant p: If we differentiate (42), we see that ∂ w l /∂t is the elliptic projection of ∂T/∂t, so Let y l = w l − T l . Subtracting (26) from (45), we get Now, for each t, we choose v = y l ðtÞ ∈ V l . Note that for any function y ∈ C 1 ð½0, T ; L 2 ðΩÞ, Thus, we get Canceling the same expression on both sides of (48), we obtain Abstract and Applied Analysis where T = kM and E is the consistency error. In this way, we prove that

Numerical Tests
We have carried out two test problems to demonstrate the performance of the given algorithm. Accuracy of the method is measured by the error norm L ∞ = kT exact − T numeric k L ∞ = kT exact − T numeric k: The equation for the demonstration is a 2D heat equation over a rectangular region: That is finding the solution of temperature in the plate as a function of time and position using the finite element method. The heat equation, which will be used for demonstration in this section, is given by where α is thermal diffusivity and f is the heat generation or source function.

Case 1: 2D Heat Equation Whose Exact Solution Is
Nonlinear. Assume that the exact solution of (67) is T x, y, t ð Þ= ty sin πx ð Þ cos π 2 y + y sin πx ð Þ cos π 2 y , 0 ≤ x, y ≤ 1, t ≥ 0: ð68Þ By using (68), we find the value of in (67). Now, the initial boundary value problem or the strong formulation is with the boundary conditions T = 0 on the boundary dΩ and initial condition Tðx, y, 0Þ = y sin ðπxÞ cos ððπ/2ÞyÞ on Ω, where f = y sin πx ð Þ 1 + απ 2 t + 1 ð Þ+ α π 2 4 t + 1 ð Þ cos π 2 y The numerical result is compared with the exact result for different values of time and number of collocation points. To demonstrate the efficiency of the method, the absolute errors are reported in some arbitrary points in Table 1. To obtain the numerical results, MATLAB software is used. The equation is solved for Δt = 0:1, α = 1, and N = 10: Its numerical and exact solution and absolute error between exact and numerical solutions can be shown in Figure 1. As seen in Table 1, the reported absolute errors are as expected, that is, of order OðΔt + l 2 Þ.

Case 2: 2D Heat Equation Whose Exact Solution Is
Linear. Assume that the exact solution of (67) is By using (72), the value of f is f = 1 in (67). Now, the initial boundary value problem or the strong formulation is with the boundary conditions Tð0, y, tÞ = y + t, Tð1, y, tÞ = 1 + y + t, Tðx, 0, tÞ = x + t, and Tðx, 1, tÞ = 1 + x + t on the boundary dΩ and initial condition Tðx, y, 0Þ = x + y on Ω, where f = 1: The numerical and exact solution and absolute error between exact and numerical solutions of example 2 can be shown in Figure 2. The numerical result is compared with the exact result for different values of time and numbers of collocation points. To demonstrate the efficiency of the used method, the absolute errors in some arbitrary points are reported in Table 2. To obtain the numerical results, MATLAB software is used. The equation is solved for Δt = 0:1, α = 1, and N = 10: Its numerical and exact solution and absolute error between exact and numerical solutions can be shown in Figure 2. As seen in Table 2, the reported absolute errors are as expected, that is, of order OðΔt + l 2 Þ.
Hence, from Tables 1 and 2, it is concluded that the absolute error between the exact solution and using FEM solution is of order OðΔt + l 2 Þ. So the finite element method is the best numerical method to solve any differential equations numerically in any type of geometries.

Conclusion
In this study, a mathematical model of a two-dimensional heat transfer equation coupled with the Darcy flow has been presented. The governing equation of a mathematical model is a system of partial differential equations and is solved using the finite element technique. After the finite element method was applied, the governing equation was then discretized into a set of ordinary differential equations. Then, backward Euler method was applied to find the numerical solution of the set of ordinary differential equations. The method was tested on two-dimensional time-dependent heat transfer in a plate, for which the exact solution is nonlinear or linear. A strong result is proven, and a numerical example 8 Abstract and Applied Analysis is provided to illustrate the convergence behavior of the problem generated by the finite element method. Further, the prior and posterior error estimates are then derived for the numerical scheme.

Data Availability
No data is available for this research.

Conflicts of Interest
The author declares that he has no conflicts of interest.