Abstract

The finite element approach was utilized in this study to solve numerically the two-dimensional time-dependent heat transfer equation coupled with the Darcy flow. 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 effectiveness of the theoretical results. The essential code development in this study was done using MATLAB computer simulation.

1. Introduction

Let be an open domain with a bounded boundary that is simply connected in , 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 and are the unknowns and 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 [28]. 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).

2. 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 , and are assigned functions and specifies a boundary partition, that is, and . is the Neumann boundary, and is the Dirichlet boundary.

An integral domain with the solution belonging to the Sobolev space and the underlying space is utilized, since the problem cannot be expressed in terms of second derivatives. The is a space that has the following definition [17]. and the Sobolev space is defined as equipped with the norm

The space for vanishing boundary values is defined as

Sobolev’s imbedding is used; for all real numbers , there exist constants and such that

Equation (10) simplifies to Poincaré’s inequality when .

2.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 with a basis function and integrating on . is set, and is sought for each such that

By putting the Darcy law, that is, , 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 with and and for all Specifying for the inner product, the variational formulation (21) and (22) can be written as

Now considering the Galerkin approximation of the variational formulation (23) and (24) for all , let be the space of vanishing Lagrange finite elements of degree on the boundary with respect to a mesh of mesh size and define and as an approximate solution, such that with and , where is an appropriate finite-dimensional space and and are convenient approximations of and in the space , respectively. The weak formulation (25) and (26) is called semidiscretization of (13) and (14).

2.3. Existence and Uniqueness

One of the objectives of the study is to construct the numerical solution of the equation of heat transfer coupled with the fluid flow equation in porous media using the finite element technique. That is the solution of (26) using the finite element technique. Since is detachable, as a result, it has a quantifiable basis Let denote the region covered by the first basis functions, The reduced problem (26) is discretized in by the square system of equations. That is, find , 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.

Theorem 1. Picard-Lindelöf Theorem: let and be separable Hilbert spaces with dense in and . Assume that there are constants with

Then, there exists a unique solution 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.

Further, let , where and Now, I am going to see if the problem meets conditions (28) and (29).

2.4. Coercivity/Positivity

We need to constrain the flow velocities to guarantee that is positive. That is what we require: and assume the bounded velocities Using the Poincaré inequality, where is some constant [19].

Now, if with , we obtain since , that is, Therefore, Hence, the positivity condition is satisfied.

2.5. Continuity/Boundedness

The second condition the boundedness follows by applying the Cauchy-Schwartz inequality: where Hence, the continuity condition is satisfied. Therefore, the solution of the weak formulation of our problem exists and is unique.

2.6. Approximate Solution

Let denote the element number in a region . To provide an algebraic interpretation of (25) and (26), a basis is introduced for , and it is observed that it suffices that (25) and (26) are verified for the basis functions in order to be satisfied by all the functions of the subspace. Moreover, as for each the solution to the finite element problem belongs to the subspace as well, thus where the and coefficients represent the unknowns of problems (25) and (26) on the element and , since the elements are triangular elements. Denoting by and the derivatives of the function and , respectively, with respect to time, the system of equations (25) and (26) on the element becomes

Then, after some implementation, systems (30) and (36) can be rewritten in the following form:

where (i) is the matrix whose elements are given by (ii) is the matrix whose elements are obtained from (19) as (iii) is the matrix whose elements are obtained from (20) as (iv) is the column vector whose elements are (v) is the column vector whose elements are

The matrix is called the local stiffness matrix, and 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 and the matrix is the global mass matrix, the matrix is the global stiffness matrix, and the vectors and are the globalized force vectors. Here, denotes the total number of nodes in the problem as a whole, and the components of and 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 .

2.7. Error Analysis
2.7.1. Prior Error Estimate

Consider the weak formulation of the governing equation (26) of heat transfer in porous 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 not directly to , but rather to an appropriate representative . For , we choose the elliptic projection of , defined by

From [20], by the finite element method for elliptic problems, we have the estimate: for some constant If we differentiate (42), we see that is the elliptic projection of , so

Now,

Let . Subtracting (26) from (45), we get

Now, for each , we choose . Note that for any function,

Thus, we get Canceling the same expression on both sides of (48), we obtain

This holds for each . Integrating over , we get

For , we have

Thus, assuming that the exact solution is sufficiently smooth and the initial data is chosen so that , we have

Combining this estimate with the elliptic estimate (43), we get an estimate on the error:

2.7.2. Posterior Error Estimate

Now, we turn to the error analysis of a fully discrete scheme that has finite elements in space and backward Euler in time. Writing for (with being the time step), the scheme is

To analyze this scheme, we proceed as we did for the semidiscrete scheme, with some extra complications coming from the time discretization. In particular, we continue to use the elliptic projection as a representative of . Thus, the consistency error is given by where the last line gives the definition of . Next, we estimate the two terms that comprise , in . First, we have by Taylor’s theorem. Next, so

Thus, we have obtained a bound on the consistency error: with

Combining with the scheme (54), we get (for )

The argument with a stability argument is concluded. Choose This becomes so and, by iteration, where and is the consistency error. In this way, we prove that

3. 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 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 is the heat generation or source function.

3.1. Case 1: 2D Heat Equation Whose Exact Solution Is Nonlinear

Assume that the exact solution of (67) is

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 on the boundary and initial condition on , where

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 , , and 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.

3.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 is in (67). Now, the initial boundary value problem or the strong formulation is with the boundary conditions and on the boundary and initial condition on , where The numerical and exact solution and absolute error between exact and numerical solutions of example 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 , , and 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.

Hence, from Tables 1 and 2, it is concluded that the absolute error between the exact solution and using FEM solution is of order . So the finite element method is the best numerical method to solve any differential equations numerically in any type of geometries.

4. 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 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.