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

Authors: Mohammed Hirpho
Journal: Abstract and Applied Analysis
 (2022)

---

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


---

## Body

## 1. 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
(1)Ss∂h∂t+∇∙q=G,(2)∂T∂t+∇T∙q=1Pe∇2T+H,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–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).


## 2. Mathematical Model Analysis

### 2.1. Initial Boundary Condition and Basic Principles

Equations (1)–(2) should be solved by specifying an initial condition:
(3)hx,y,0=h0x,y,x,y∈Ω,Tx,y,0=T0x,y,x,y∈Ω,as well as boundary conditions, which could also take the following form:
(4)Tx,y,t=px,y,t,hx,y,t=qx,y,t,x,y∈∂ΩD,t>0,∂T∂n=rx,y,t,∂h∂n=sx,y,t,x,y∈∂ΩD,t>0,where T0,h0,p,q,r, and s are assigned functions and ∂ΩD,∂ΩN 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 spaceH1Ω and the underlying space L2Ω is utilized, since the problem cannot be expressed in terms of second derivatives. The L2Ω is a space that has the following definition [17].
(5)L2Ω=v:∫Ω.vx,y2dΩ<∞,and the Sobolev space H1Ω is defined as
(6)H1Ω=v∈L2Ω:∂v∂x,∂v∂y∈L2Ω,equipped with the norm
(7)vH1Ω=∑vH1Ω21/2.The space for vanishing boundary values is defined as(8)H01Ω=v∈H1Ω:v∂ΩD=px,y.Sobolev’s imbedding is used; for all real numbersp≥1, there exist constants Sp0 and Sp such that
(9)∀v∈H1Ω,vLpΩ≥SpvH1Ω,(10)∀v∈H01Ω,vLpΩ≥Sp0vH1Ω.Equation (10) simplifies to Poincaré’s inequality when p=2.


### 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 t>0 with a basis function v=vx,y and integrating on Ω. V=H1Ω is set, and ht,Tt∈V is sought for each t>0 such that
(11)∫Ω.Ss∂h∂tvdΩ+∇∙qvdΩ=∫Ω.GvdΩ,(12)∫Ω.∂T∂tvdΩ+∫Ω.∇T∙qvdΩ=1Pe∫Ω.∇2TvdΩ+∫Ω.HvdΩ.By putting the Darcy law, that is,q=−K∇h, into (11) and (12), we obtain
(13)∫Ω.Ss∂h∂tvdΩ−∫Ω.K∇2hvdΩ=∫Ω.GvdΩ,(14)∫Ω.∂T∂tvdΩ−∫Ω.∇T∙K∇h+1Pe∇2TvdΩ=∫Ω.HvdΩ.Then, by using Green’s theorem stated in [18], equations (13) and (14) can be
(15)∫Ω.Ss∂h∂tvdΩ+∫Ω.K∇h∇vdΩ−∫∂Ωn.v∂h∂nds=∫Ω.GvdΩ,(16)∫Ω.∂T∂tvdΩ−∫Ω.∇T∙K∇hvdΩ+1Pe∫Ω.∇T∇vdΩ−∫∂Ωn.v∂T∂nds=∫Ω.HvdΩ.In (16),
(17)−∫Ω∇T·K∇hvdΩ=∫Ω∇T·qvdΩ,q=−K∇h,⇒−∫Ω∇T·K∇hvdΩ=∫Ω∇T·qvdΩ=∫Ωq∇Tv+∇qTvdΩ,since we considered basis functions to be a linear basis function (e.g., ∇q=∇−K∇h=0),
(18)∫Ωq∇Tv+∇qTvdΩ=∫Ωq∇TvdΩ.Assume(19)aht,v=∫Ω.K∇h∇vdΩ,(20)bTt,v=−∫Ω.∇T∙K∇hvdΩ+1Pe∫Ω.∇T∇vdΩ.Bringing (19) and (20) into (15) and (16), respectively, we obtain
(21)∫Ω.Ss∂h∂tvdΩ+aht,v=∫Ω.GvdΩ+∫∂Ωnv∂h∂nds,(22)∫Ω.∂T∂tvdΩ+bTt,v=∫Ω.HvdΩ+∫∂Ωnv∂T∂nds,with T0=T0 and h0=h0 and for all v∈V. Specifying .,. for the L2Ω inner product, the variational formulation (21) and (22) can be written as
(23)∂h∂t,v+aht,v=G,v+∫∂Ωneϕisx,yds,v∈V,(24)∂T∂t,v+bTt,v=H,v+∫∂Ωneϕirx,yds,v∈V.Now considering the Galerkin approximation of the variational formulation (23) and (24) for all t>0, let Vl 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 hl and Tl as an approximate solution, hlt,Tlt∈Vl such that
(25)∂hl∂t,vl+ahlt,vl=G,vl+∫∂Ωneϕisx,yds,vl∈Vl,(26)∂Tl∂t,vl+bTlt,vl=H,vl+∫∂Ωneϕirx,yds,vl∈Vl,with Tl0=T0l and hl0=h0l, where Vl⊂V is an appropriate finite-dimensional space and T0l and h0l are convenient approximations of T0 and h0 in the space Vl, 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 H1Ω is detachable, as a result, it has a quantifiable basis ϕii≥1. Let Vl denote the region covered by the first l basis functions, ϕi1≤i≤l. The reduced problem (26) is discretized in Vm by the square system of equations. That is, find Tl=∑1≤i≤lTiϕi∈Vl, the solution of
(27)∂Tlt∂t,vl+bTlt,vl=H,vl,Tl0,vl=T0,vl,vl∈Vl.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: letV and H be separable Hilbert spaces with V⊂H,V dense in H and ∥v∥≤c∥v∥V. Assume that there are constants μ2≥μ1>0 with
(28)bv,v≥μ1v2V,∀v∈V,(29)bw,v≤μ2wVvV,∀w,v∈V.Then, there exists a unique solutionT∈L20,Tf,V of the initial value problem:
(30)∂Tt∂t,v+bTt,v=H,v,T0=T0,v∈V.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, letbw,v=cw,v+dw,v, where cw,v=1/Pe∫Ω∇w∇vdΩ and dw,v=∫Ωq∇wvdΩ. 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 thatbv,v is positive. That is what we require:
(31)dv,v≤bv,v,and assume the bounded velocities q∞. Using the Poincaré inequality, vH01≤CΩvH1, where CΩ is some constant [19].Now, ifCΩq∞≤D1/Pe, with D<1, we obtain
(32)bv,v=1Pe∫Ω∇v∇vdΩ+∫Ωq∇vvdΩ≥1PevH1vH1−q∞vH1vH01,≥1Pe−CΩq∞vH1vH1≥1Pe1−DvH12,since D<1, that is, μ1=1/Pe1−D>0. Therefore, bv,v≥μ1vH12. Hence, the positivity condition is satisfied.


### 2.5. Continuity/Boundedness

The second condition the boundedness follows by applying the Cauchy-Schwartz inequality:(33)bw,v=1Pe∫Ω∇w∇vdΩ+∫Ωq∇wvdΩ≤1PewH1vH1+q∞wH1vH01,≤1Pe+CΩq∞wH1vH1≤μ2wH1vH1,where μ2=1/Pe+CΩq∞>0. Hence, the continuity condition is satisfied. Therefore, the solution of the weak formulation of our problem exists and is unique.


### 2.6. Approximate Solution

Lete denote the element number in a region Ω. To provide an algebraic interpretation of (25) and (26), a basis ϕj is introduced for Vl, 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 t>0 the solution to the finite element problem belongs to the subspace as well, thus
(34)Tlex,y,t=∑j=1NlTjetϕjx,y,hlex,y,t=∑j=1Nlhjetϕjx,y,where the Tjet and hjet coefficients represent the unknowns of problems (25) and (26) on the element e and Nl=3, since the elements are triangular elements. Denoting by Tėjt and hėjt the derivatives of the function Tjet and hjet, respectively, with respect to time, the system of equations (25) and (26) on the element e becomes
(35)∫Ωe∑j=1NlhėjtϕjϕidΩe+a∑j=1Nlhjetϕj,ϕi=∫ΩeGϕidΩe+∫∂Ωneϕisx,yds,i=1,2,⋯,Nl,(36)∫Ωe∑j=1NlTėjtSsϕjϕidΩe+b∑j=1NlTjetϕj,ϕi=∫ΩeHϕidΩe+∫∂Ωneϕirx,yds,i=1,2,⋯,Nl.Then, after some implementation, systems (30) and (36) can be rewritten in the following form:
(37)SsMeḣe+Aehe=ge,(38)MeṪe+BeTe=re,where(i)
Me is the 3×3 matrix whose elements are given by mjie=∫ΩeϕjϕidΩe(ii)
Ae is the 3×3 matrix whose elements are obtained from (19) as ajie=∫ΩeK∇ϕj∇ϕidΩe(iii)
Be is the 3×3 matrix whose elements are obtained from (20) as bjie=∫Ωeq∇ϕjϕidΩe+1/Pe∫Ωe∇ϕj∇ϕidΩe(iv)
ge is the column vector whose elements are gie=∫ΩeGϕidΩe+∫∂Ωneϕisx,yds(v)
re is the column vector whose elements are rie=∫ΩeHϕidΩe+∫∂Ωneϕirx,ydsThe matrixAe is called the local stiffness matrix, and Me 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:
(39)SsMḣ+Ah=g,(40)MṪ+BT=r,in which h=h1,h2,⋯,hNT,T=T1,T2,⋯,TNT 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:
(41)SsMhk+1−hkΔt+Chk+1=gk+1,MTk+1−TkΔt+CTk+1=rk+1,which is first order accurate with respect to Δt=tk+1−tk.


### 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 Tl not directly to T, but rather to an appropriate representative wl∈C10,Tf,Vl. For wl, we choose the elliptic projection of T, defined by
(42)bwl,v=bT,v,v∈Vl,0≤t≤Tf.From [20], by the finite element method for elliptic problems, we have the L2 estimate:
(43)Tt−wlt≤plr+1Ttr+1,0≤t≤Tf,for some constant p. If we differentiate (42), we see that ∂wl/∂t is the elliptic projection of ∂T/∂t, so
(44)∂T∂tt−∂wl∂tt≤plr+1∂T∂ttr+1,0≤t≤Tf.Now,(45)∂wl∂t,v+bwl,v=∂wl∂t,v+bT,v=∂wl−T∂t,v+H,v,v∈Vl,0≤t≤TfLetyl=wl−Tl. Subtracting (26) from (45), we get
(46)∂yl∂t,v+byl,v=∂wl−T∂t,v,v∈Vl,0≤t≤Tf.Now, for eacht, we choose v=ylt∈Vl. Note that for any functiony∈C10,T;L2Ω,
(47)∥y∥ddt∥y∥=12ddty2=∂y∂t,y.Thus, we get(48)∥yl∥ddt∥yl∥+byl,yl=∂wl−T∂t,yl≤∂wl−T∂t∥yl∥.Canceling the same expression on both sides of (48), we obtain
(49)ddtyl≤∂wl−T∂t≤plr+1∂T∂ttr+1.This holds for eacht. Integrating over 0,t, we get
(50)ylt≤wlt−Tt≤yl0+plr+1TtL10,T;Hr+1Ω.Foryl0, we have
(51)yl0=wl0−Tl0≤wl0−T0+T0−Tl0≤plr+1T0r+1∥+∥T0−Tl0∥.Thus, assuming that the exact solution is sufficiently smooth and the initial dataTl0 is chosen so that ∥T0−Tl0∥=Olr+1, we have
(52)ylL∞0,T;L2Ω=Olr+1.Combining this estimate with the elliptic estimate (43), we get an estimate on the error:
(53)T−TlL∞0,T;L2Ω=Olr+1.


#### 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. WritingTlj for Tl·,jk (with k being the time step), the scheme is
(54)Tlj+1−Tljk,v+bTlj+1,v=Hj+1,v,v∈Vl,j=0,1,2,⋯,M.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 projectionwl as a representative of T. Thus, the consistency error is given by
(55)E=wlj+1−wljk,v+bwlj+1,v−Hj+1,v=Tj+1−Tjk,v+bTj+1,v−Hj+1,v+wlj+1−Tj+1−wlj−Tjk,v=Tj+1−Tjk−∂Tj+1∂t,v+wlj+1−Tj+1−wlj−Tjk,v=zj,v,where the last line gives the definition of zj. Next, we estimate the two terms that comprise zj, in L2. First, we have
(56)Tj+1−Tjk−∂Tj+1∂t≤k2∂2T∂t2L∞L2,by Taylor’s theorem. Next,
(57)wlj+1−Tj+1−wlj−Tjk=1k∫jkj+1k∂∂twls−Tsds,so
(58)wlj+1−Tj+1−wlj−Tjk=plr+1∂T∂tL∞jk,j+k,Hr+1Ω.Thus, we have obtained a bound on the consistency error:(59)E=wlj+1−wljk,v+bwlj+1,v−Hj+1,v=zj,v,v∈Vl,j=0,1,2,⋯,M,with
(60)zj=pk∂2T∂t2L∞0,Tf,Hr+1Ω+lr+1∂T∂tL∞0,Tf,Hr+1Ω.Combining with the scheme (54), we get (for yl=wl−Tl)
(61)ylj+1−yljk,v+bylj+1,v=zj,v,v∈Vl.The argument with a stability argument is concluded. Choosev=ylj+1∈Vl. This becomes
(62)ylj+12+kbylj+1,ylj+1=ylj+kzj,ylj+1,so
(63)ylj+1≤ylj+kE,j=0,1,⋯,M,and, by iteration,
(64)max0≤j≤Mylj≤yl0+TE,j=0,1,⋯,M,where T=kM and E is the consistency error. In this way, we prove that
(65)max0≤j≤MTj−Tlj=Ok+lr+1,j=0,1,⋯,M.


## 2.1. Initial Boundary Condition and Basic Principles

Equations (1)–(2) should be solved by specifying an initial condition:
(3)hx,y,0=h0x,y,x,y∈Ω,Tx,y,0=T0x,y,x,y∈Ω,as well as boundary conditions, which could also take the following form:
(4)Tx,y,t=px,y,t,hx,y,t=qx,y,t,x,y∈∂ΩD,t>0,∂T∂n=rx,y,t,∂h∂n=sx,y,t,x,y∈∂ΩD,t>0,where T0,h0,p,q,r, and s are assigned functions and ∂ΩD,∂ΩN 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 spaceH1Ω and the underlying space L2Ω is utilized, since the problem cannot be expressed in terms of second derivatives. The L2Ω is a space that has the following definition [17].
(5)L2Ω=v:∫Ω.vx,y2dΩ<∞,and the Sobolev space H1Ω is defined as
(6)H1Ω=v∈L2Ω:∂v∂x,∂v∂y∈L2Ω,equipped with the norm
(7)vH1Ω=∑vH1Ω21/2.The space for vanishing boundary values is defined as(8)H01Ω=v∈H1Ω:v∂ΩD=px,y.Sobolev’s imbedding is used; for all real numbersp≥1, there exist constants Sp0 and Sp such that
(9)∀v∈H1Ω,vLpΩ≥SpvH1Ω,(10)∀v∈H01Ω,vLpΩ≥Sp0vH1Ω.Equation (10) simplifies to Poincaré’s inequality when p=2.


## 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 t>0 with a basis function v=vx,y and integrating on Ω. V=H1Ω is set, and ht,Tt∈V is sought for each t>0 such that
(11)∫Ω.Ss∂h∂tvdΩ+∇∙qvdΩ=∫Ω.GvdΩ,(12)∫Ω.∂T∂tvdΩ+∫Ω.∇T∙qvdΩ=1Pe∫Ω.∇2TvdΩ+∫Ω.HvdΩ.By putting the Darcy law, that is,q=−K∇h, into (11) and (12), we obtain
(13)∫Ω.Ss∂h∂tvdΩ−∫Ω.K∇2hvdΩ=∫Ω.GvdΩ,(14)∫Ω.∂T∂tvdΩ−∫Ω.∇T∙K∇h+1Pe∇2TvdΩ=∫Ω.HvdΩ.Then, by using Green’s theorem stated in [18], equations (13) and (14) can be
(15)∫Ω.Ss∂h∂tvdΩ+∫Ω.K∇h∇vdΩ−∫∂Ωn.v∂h∂nds=∫Ω.GvdΩ,(16)∫Ω.∂T∂tvdΩ−∫Ω.∇T∙K∇hvdΩ+1Pe∫Ω.∇T∇vdΩ−∫∂Ωn.v∂T∂nds=∫Ω.HvdΩ.In (16),
(17)−∫Ω∇T·K∇hvdΩ=∫Ω∇T·qvdΩ,q=−K∇h,⇒−∫Ω∇T·K∇hvdΩ=∫Ω∇T·qvdΩ=∫Ωq∇Tv+∇qTvdΩ,since we considered basis functions to be a linear basis function (e.g., ∇q=∇−K∇h=0),
(18)∫Ωq∇Tv+∇qTvdΩ=∫Ωq∇TvdΩ.Assume(19)aht,v=∫Ω.K∇h∇vdΩ,(20)bTt,v=−∫Ω.∇T∙K∇hvdΩ+1Pe∫Ω.∇T∇vdΩ.Bringing (19) and (20) into (15) and (16), respectively, we obtain
(21)∫Ω.Ss∂h∂tvdΩ+aht,v=∫Ω.GvdΩ+∫∂Ωnv∂h∂nds,(22)∫Ω.∂T∂tvdΩ+bTt,v=∫Ω.HvdΩ+∫∂Ωnv∂T∂nds,with T0=T0 and h0=h0 and for all v∈V. Specifying .,. for the L2Ω inner product, the variational formulation (21) and (22) can be written as
(23)∂h∂t,v+aht,v=G,v+∫∂Ωneϕisx,yds,v∈V,(24)∂T∂t,v+bTt,v=H,v+∫∂Ωneϕirx,yds,v∈V.Now considering the Galerkin approximation of the variational formulation (23) and (24) for all t>0, let Vl 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 hl and Tl as an approximate solution, hlt,Tlt∈Vl such that
(25)∂hl∂t,vl+ahlt,vl=G,vl+∫∂Ωneϕisx,yds,vl∈Vl,(26)∂Tl∂t,vl+bTlt,vl=H,vl+∫∂Ωneϕirx,yds,vl∈Vl,with Tl0=T0l and hl0=h0l, where Vl⊂V is an appropriate finite-dimensional space and T0l and h0l are convenient approximations of T0 and h0 in the space Vl, 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 H1Ω is detachable, as a result, it has a quantifiable basis ϕii≥1. Let Vl denote the region covered by the first l basis functions, ϕi1≤i≤l. The reduced problem (26) is discretized in Vm by the square system of equations. That is, find Tl=∑1≤i≤lTiϕi∈Vl, the solution of
(27)∂Tlt∂t,vl+bTlt,vl=H,vl,Tl0,vl=T0,vl,vl∈Vl.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: letV and H be separable Hilbert spaces with V⊂H,V dense in H and ∥v∥≤c∥v∥V. Assume that there are constants μ2≥μ1>0 with
(28)bv,v≥μ1v2V,∀v∈V,(29)bw,v≤μ2wVvV,∀w,v∈V.Then, there exists a unique solutionT∈L20,Tf,V of the initial value problem:
(30)∂Tt∂t,v+bTt,v=H,v,T0=T0,v∈V.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, letbw,v=cw,v+dw,v, where cw,v=1/Pe∫Ω∇w∇vdΩ and dw,v=∫Ωq∇wvdΩ. 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 thatbv,v is positive. That is what we require:
(31)dv,v≤bv,v,and assume the bounded velocities q∞. Using the Poincaré inequality, vH01≤CΩvH1, where CΩ is some constant [19].Now, ifCΩq∞≤D1/Pe, with D<1, we obtain
(32)bv,v=1Pe∫Ω∇v∇vdΩ+∫Ωq∇vvdΩ≥1PevH1vH1−q∞vH1vH01,≥1Pe−CΩq∞vH1vH1≥1Pe1−DvH12,since D<1, that is, μ1=1/Pe1−D>0. Therefore, bv,v≥μ1vH12. Hence, the positivity condition is satisfied.


## 2.5. Continuity/Boundedness

The second condition the boundedness follows by applying the Cauchy-Schwartz inequality:(33)bw,v=1Pe∫Ω∇w∇vdΩ+∫Ωq∇wvdΩ≤1PewH1vH1+q∞wH1vH01,≤1Pe+CΩq∞wH1vH1≤μ2wH1vH1,where μ2=1/Pe+CΩq∞>0. Hence, the continuity condition is satisfied. Therefore, the solution of the weak formulation of our problem exists and is unique.


## 2.6. Approximate Solution

Lete denote the element number in a region Ω. To provide an algebraic interpretation of (25) and (26), a basis ϕj is introduced for Vl, 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 t>0 the solution to the finite element problem belongs to the subspace as well, thus
(34)Tlex,y,t=∑j=1NlTjetϕjx,y,hlex,y,t=∑j=1Nlhjetϕjx,y,where the Tjet and hjet coefficients represent the unknowns of problems (25) and (26) on the element e and Nl=3, since the elements are triangular elements. Denoting by Tėjt and hėjt the derivatives of the function Tjet and hjet, respectively, with respect to time, the system of equations (25) and (26) on the element e becomes
(35)∫Ωe∑j=1NlhėjtϕjϕidΩe+a∑j=1Nlhjetϕj,ϕi=∫ΩeGϕidΩe+∫∂Ωneϕisx,yds,i=1,2,⋯,Nl,(36)∫Ωe∑j=1NlTėjtSsϕjϕidΩe+b∑j=1NlTjetϕj,ϕi=∫ΩeHϕidΩe+∫∂Ωneϕirx,yds,i=1,2,⋯,Nl.Then, after some implementation, systems (30) and (36) can be rewritten in the following form:
(37)SsMeḣe+Aehe=ge,(38)MeṪe+BeTe=re,where(i)
Me is the 3×3 matrix whose elements are given by mjie=∫ΩeϕjϕidΩe(ii)
Ae is the 3×3 matrix whose elements are obtained from (19) as ajie=∫ΩeK∇ϕj∇ϕidΩe(iii)
Be is the 3×3 matrix whose elements are obtained from (20) as bjie=∫Ωeq∇ϕjϕidΩe+1/Pe∫Ωe∇ϕj∇ϕidΩe(iv)
ge is the column vector whose elements are gie=∫ΩeGϕidΩe+∫∂Ωneϕisx,yds(v)
re is the column vector whose elements are rie=∫ΩeHϕidΩe+∫∂Ωneϕirx,ydsThe matrixAe is called the local stiffness matrix, and Me 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:
(39)SsMḣ+Ah=g,(40)MṪ+BT=r,in which h=h1,h2,⋯,hNT,T=T1,T2,⋯,TNT 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:
(41)SsMhk+1−hkΔt+Chk+1=gk+1,MTk+1−TkΔt+CTk+1=rk+1,which is first order accurate with respect to Δt=tk+1−tk.


## 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 Tl not directly to T, but rather to an appropriate representative wl∈C10,Tf,Vl. For wl, we choose the elliptic projection of T, defined by
(42)bwl,v=bT,v,v∈Vl,0≤t≤Tf.From [20], by the finite element method for elliptic problems, we have the L2 estimate:
(43)Tt−wlt≤plr+1Ttr+1,0≤t≤Tf,for some constant p. If we differentiate (42), we see that ∂wl/∂t is the elliptic projection of ∂T/∂t, so
(44)∂T∂tt−∂wl∂tt≤plr+1∂T∂ttr+1,0≤t≤Tf.Now,(45)∂wl∂t,v+bwl,v=∂wl∂t,v+bT,v=∂wl−T∂t,v+H,v,v∈Vl,0≤t≤TfLetyl=wl−Tl. Subtracting (26) from (45), we get
(46)∂yl∂t,v+byl,v=∂wl−T∂t,v,v∈Vl,0≤t≤Tf.Now, for eacht, we choose v=ylt∈Vl. Note that for any functiony∈C10,T;L2Ω,
(47)∥y∥ddt∥y∥=12ddty2=∂y∂t,y.Thus, we get(48)∥yl∥ddt∥yl∥+byl,yl=∂wl−T∂t,yl≤∂wl−T∂t∥yl∥.Canceling the same expression on both sides of (48), we obtain
(49)ddtyl≤∂wl−T∂t≤plr+1∂T∂ttr+1.This holds for eacht. Integrating over 0,t, we get
(50)ylt≤wlt−Tt≤yl0+plr+1TtL10,T;Hr+1Ω.Foryl0, we have
(51)yl0=wl0−Tl0≤wl0−T0+T0−Tl0≤plr+1T0r+1∥+∥T0−Tl0∥.Thus, assuming that the exact solution is sufficiently smooth and the initial dataTl0 is chosen so that ∥T0−Tl0∥=Olr+1, we have
(52)ylL∞0,T;L2Ω=Olr+1.Combining this estimate with the elliptic estimate (43), we get an estimate on the error:
(53)T−TlL∞0,T;L2Ω=Olr+1.


### 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. WritingTlj for Tl·,jk (with k being the time step), the scheme is
(54)Tlj+1−Tljk,v+bTlj+1,v=Hj+1,v,v∈Vl,j=0,1,2,⋯,M.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 projectionwl as a representative of T. Thus, the consistency error is given by
(55)E=wlj+1−wljk,v+bwlj+1,v−Hj+1,v=Tj+1−Tjk,v+bTj+1,v−Hj+1,v+wlj+1−Tj+1−wlj−Tjk,v=Tj+1−Tjk−∂Tj+1∂t,v+wlj+1−Tj+1−wlj−Tjk,v=zj,v,where the last line gives the definition of zj. Next, we estimate the two terms that comprise zj, in L2. First, we have
(56)Tj+1−Tjk−∂Tj+1∂t≤k2∂2T∂t2L∞L2,by Taylor’s theorem. Next,
(57)wlj+1−Tj+1−wlj−Tjk=1k∫jkj+1k∂∂twls−Tsds,so
(58)wlj+1−Tj+1−wlj−Tjk=plr+1∂T∂tL∞jk,j+k,Hr+1Ω.Thus, we have obtained a bound on the consistency error:(59)E=wlj+1−wljk,v+bwlj+1,v−Hj+1,v=zj,v,v∈Vl,j=0,1,2,⋯,M,with
(60)zj=pk∂2T∂t2L∞0,Tf,Hr+1Ω+lr+1∂T∂tL∞0,Tf,Hr+1Ω.Combining with the scheme (54), we get (for yl=wl−Tl)
(61)ylj+1−yljk,v+bylj+1,v=zj,v,v∈Vl.The argument with a stability argument is concluded. Choosev=ylj+1∈Vl. This becomes
(62)ylj+12+kbylj+1,ylj+1=ylj+kzj,ylj+1,so
(63)ylj+1≤ylj+kE,j=0,1,⋯,M,and, by iteration,
(64)max0≤j≤Mylj≤yl0+TE,j=0,1,⋯,M,where T=kM and E is the consistency error. In this way, we prove that
(65)max0≤j≤MTj−Tlj=Ok+lr+1,j=0,1,⋯,M.


## 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 Tl not directly to T, but rather to an appropriate representative wl∈C10,Tf,Vl. For wl, we choose the elliptic projection of T, defined by
(42)bwl,v=bT,v,v∈Vl,0≤t≤Tf.From [20], by the finite element method for elliptic problems, we have the L2 estimate:
(43)Tt−wlt≤plr+1Ttr+1,0≤t≤Tf,for some constant p. If we differentiate (42), we see that ∂wl/∂t is the elliptic projection of ∂T/∂t, so
(44)∂T∂tt−∂wl∂tt≤plr+1∂T∂ttr+1,0≤t≤Tf.Now,(45)∂wl∂t,v+bwl,v=∂wl∂t,v+bT,v=∂wl−T∂t,v+H,v,v∈Vl,0≤t≤TfLetyl=wl−Tl. Subtracting (26) from (45), we get
(46)∂yl∂t,v+byl,v=∂wl−T∂t,v,v∈Vl,0≤t≤Tf.Now, for eacht, we choose v=ylt∈Vl. Note that for any functiony∈C10,T;L2Ω,
(47)∥y∥ddt∥y∥=12ddty2=∂y∂t,y.Thus, we get(48)∥yl∥ddt∥yl∥+byl,yl=∂wl−T∂t,yl≤∂wl−T∂t∥yl∥.Canceling the same expression on both sides of (48), we obtain
(49)ddtyl≤∂wl−T∂t≤plr+1∂T∂ttr+1.This holds for eacht. Integrating over 0,t, we get
(50)ylt≤wlt−Tt≤yl0+plr+1TtL10,T;Hr+1Ω.Foryl0, we have
(51)yl0=wl0−Tl0≤wl0−T0+T0−Tl0≤plr+1T0r+1∥+∥T0−Tl0∥.Thus, assuming that the exact solution is sufficiently smooth and the initial dataTl0 is chosen so that ∥T0−Tl0∥=Olr+1, we have
(52)ylL∞0,T;L2Ω=Olr+1.Combining this estimate with the elliptic estimate (43), we get an estimate on the error:
(53)T−TlL∞0,T;L2Ω=Olr+1.


## 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. WritingTlj for Tl·,jk (with k being the time step), the scheme is
(54)Tlj+1−Tljk,v+bTlj+1,v=Hj+1,v,v∈Vl,j=0,1,2,⋯,M.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 projectionwl as a representative of T. Thus, the consistency error is given by
(55)E=wlj+1−wljk,v+bwlj+1,v−Hj+1,v=Tj+1−Tjk,v+bTj+1,v−Hj+1,v+wlj+1−Tj+1−wlj−Tjk,v=Tj+1−Tjk−∂Tj+1∂t,v+wlj+1−Tj+1−wlj−Tjk,v=zj,v,where the last line gives the definition of zj. Next, we estimate the two terms that comprise zj, in L2. First, we have
(56)Tj+1−Tjk−∂Tj+1∂t≤k2∂2T∂t2L∞L2,by Taylor’s theorem. Next,
(57)wlj+1−Tj+1−wlj−Tjk=1k∫jkj+1k∂∂twls−Tsds,so
(58)wlj+1−Tj+1−wlj−Tjk=plr+1∂T∂tL∞jk,j+k,Hr+1Ω.Thus, we have obtained a bound on the consistency error:(59)E=wlj+1−wljk,v+bwlj+1,v−Hj+1,v=zj,v,v∈Vl,j=0,1,2,⋯,M,with
(60)zj=pk∂2T∂t2L∞0,Tf,Hr+1Ω+lr+1∂T∂tL∞0,Tf,Hr+1Ω.Combining with the scheme (54), we get (for yl=wl−Tl)
(61)ylj+1−yljk,v+bylj+1,v=zj,v,v∈Vl.The argument with a stability argument is concluded. Choosev=ylj+1∈Vl. This becomes
(62)ylj+12+kbylj+1,ylj+1=ylj+kzj,ylj+1,so
(63)ylj+1≤ylj+kE,j=0,1,⋯,M,and, by iteration,
(64)max0≤j≤Mylj≤yl0+TE,j=0,1,⋯,M,where T=kM and E is the consistency error. In this way, we prove that
(65)max0≤j≤MTj−Tlj=Ok+lr+1,j=0,1,⋯,M.


## 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 normL∞=Texact−TnumericL∞=Texact−Tnumeric. The equation for the demonstration is a 2D heat equation over a rectangular region:
(66)Ω=x,y:0≤x≤1,0≤y≤1.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(67)∂T∂t=α∂2T∂x2+∂2T∂y2+f,0≤x,y≤1,t≥0,where α is thermal diffusivity and f 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
(68)Tx,y,t=tysinπxcosπ2y+ysinπxcosπ2y,0≤x,y≤1,t≥0.By using (68), we find the value of
(69)f=ysinπx1+απ2t+1+απ24t+1cosπ2y+απt+1sinπ2y,in (67). Now, the initial boundary value problem or the strong formulation is
(70)∂T∂t=α∂2T∂x2+∂2T∂y2+f,0≤x,y≤1,t≥0,with the boundary conditions T=0 on the boundary dΩ and initial condition Tx,y,0=ysinπxcosπ/2y on Ω, where
(71)f=ysinπx1+απ2t+1+απ24t+1cosπ2y+απt+1sinπ2y.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 Table1. 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 orderOΔt+l2.Table 1
Comparison of the exact and the numerical solutions using the finite element method forΔt=0.1, α=1, and N=10.
txyExact solutionFEM solutionAbsolute error0000010.30.30.41090.41430.00340.60.60.63730.64980.01251100000000150.30.33.13573.57540.43980.60.64.86345.62840.76511100000000500.30.39.947610.07800.13040.60.615.428915.95240.523511000000001000.30.319.678918.88040.79850.60.630.522330.02400.498411000Figure 1
The numerical solution, exact solution, and absolute error between them, respectively, when the exact solution is nonlinear.
(a)
Numerical solution(b)
Exact solution(c)
Absolute error


### 3.2. Case 2: 2D Heat Equation Whose Exact Solution Is Linear

Assume that the exact solution of (67) is
(72)Tx,y,t=x+y+t,0≤x,y≤1,t≥0.By using (72), the value of f is f=1 in (67). Now, the initial boundary value problem or the strong formulation is
(73)∂T∂t=α∂2T∂x2+∂2T∂y2+f,0≤x,y≤1,t≥0,with the boundary conditions T0,y,t=y+t,T1,y,t=1+y+t,Tx,0,t=x+t, and Tx,1,t=1+x+t on the boundary dΩ and initial condition Tx,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.Figure 2
The numerical solution, exact solution, and absolute error between them, respectively, when the exact solution is linear.
(a)
Numerical solution(b)
Exact solution(c)
Absolute errorThe 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 Table2. 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 orderOΔt+l2.Table 2
Comparison of the exact and the numerical solutions using the finite element method forΔt=0.1, α=1, and N=10.
txyFE solutionExact solutionAbsolute error002.00002.0000010.30.32.05442.00000.05440.60.62.06822.00000.0682112.00002.000000016.000016.00000150.30.316.054416.00000.05440.60.616.068216.00000.06821116.000016.000000051.000051.00000500.30.351.054451.00000.05440.60.651.068251.00000.06821151.000051.0000000101.0000101.000001000.30.3101.0544101.00000.05440.60.6101.0682101.00000.068211101.0000101.00000Hence, from Tables1 and 2, it is concluded that the absolute error between the exact solution and using FEM solution is of order OΔt+l2. So the finite element method is the best numerical method to solve any differential equations numerically in any type of geometries.


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

Assume that the exact solution of (67) is
(68)Tx,y,t=tysinπxcosπ2y+ysinπxcosπ2y,0≤x,y≤1,t≥0.By using (68), we find the value of
(69)f=ysinπx1+απ2t+1+απ24t+1cosπ2y+απt+1sinπ2y,in (67). Now, the initial boundary value problem or the strong formulation is
(70)∂T∂t=α∂2T∂x2+∂2T∂y2+f,0≤x,y≤1,t≥0,with the boundary conditions T=0 on the boundary dΩ and initial condition Tx,y,0=ysinπxcosπ/2y on Ω, where
(71)f=ysinπx1+απ2t+1+απ24t+1cosπ2y+απt+1sinπ2y.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 Table1. 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 orderOΔt+l2.Table 1
Comparison of the exact and the numerical solutions using the finite element method forΔt=0.1, α=1, and N=10.
txyExact solutionFEM solutionAbsolute error0000010.30.30.41090.41430.00340.60.60.63730.64980.01251100000000150.30.33.13573.57540.43980.60.64.86345.62840.76511100000000500.30.39.947610.07800.13040.60.615.428915.95240.523511000000001000.30.319.678918.88040.79850.60.630.522330.02400.498411000Figure 1
The numerical solution, exact solution, and absolute error between them, respectively, when the exact solution is nonlinear.
(a)
Numerical solution(b)
Exact solution(c)
Absolute error


## 3.2. Case 2: 2D Heat Equation Whose Exact Solution Is Linear

Assume that the exact solution of (67) is
(72)Tx,y,t=x+y+t,0≤x,y≤1,t≥0.By using (72), the value of f is f=1 in (67). Now, the initial boundary value problem or the strong formulation is
(73)∂T∂t=α∂2T∂x2+∂2T∂y2+f,0≤x,y≤1,t≥0,with the boundary conditions T0,y,t=y+t,T1,y,t=1+y+t,Tx,0,t=x+t, and Tx,1,t=1+x+t on the boundary dΩ and initial condition Tx,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.Figure 2
The numerical solution, exact solution, and absolute error between them, respectively, when the exact solution is linear.
(a)
Numerical solution(b)
Exact solution(c)
Absolute errorThe 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 Table2. 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 orderOΔt+l2.Table 2
Comparison of the exact and the numerical solutions using the finite element method forΔt=0.1, α=1, and N=10.
txyFE solutionExact solutionAbsolute error002.00002.0000010.30.32.05442.00000.05440.60.62.06822.00000.0682112.00002.000000016.000016.00000150.30.316.054416.00000.05440.60.616.068216.00000.06821116.000016.000000051.000051.00000500.30.351.054451.00000.05440.60.651.068251.00000.06821151.000051.0000000101.0000101.000001000.30.3101.0544101.00000.05440.60.6101.0682101.00000.068211101.0000101.00000Hence, from Tables1 and 2, it is concluded that the absolute error between the exact solution and using FEM solution is of order OΔt+l2. 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.



---

*Source: 5108445-2022-04-05.xml*