# A Leap-Frog Finite Difference Method for Strongly Coupled System from Sweat Transport in Porous Textile Media

Authors: Qian Zhang; Chao Huang
Journal: Journal of Mathematics
 (2019)

---

## Abstract

In this paper, we present an uncoupled leap-frog finite difference method for the system of equations arising from sweat transport through porous textile media. Based on physical mechanisms, the sweat transport can be viewed as the multicomponent flow that coupled the heat and moisture transfer, such that the system is nonlinear and strongly coupled. The leap-frog method is proposed to solve this system, with the second order accuracy in both spatial and temporal directions. We prove the existence and uniqueness of the solution to the system with optimal error estimates in the discreteL2 norm. Numerical simulations are presented and analyzed, respectively.


---

## Body

## 1. Introduction

Single/multicomponent flow in porous textile media attracted considerable attention in the last several decades. See [1–4] for the single-component models and [5–9] for the multicomponent models. In this paper, we study the multicomponent sweat transport coupled with vapor and heat in porous textile media. In [10], Ye et al. proposed a quasi-steady-state single-component model which consists of a steady-state air equation and dynamic state equations for other components. Under certain conditions, the multicomponent model reduces to a new single-component model, and the physical process can be viewed as sweat transport (vapor and heat flow) governed by the conservation of mass and energy: (1)∂∂tϵC+∂∂xugϵC=-Γce,(2)∂∂tCvtT+∂∂xϵugCmgCT-∂∂xκ∂T∂x=λMwΓce,where ϵ is the porosity of the media, C is the vapor concentration, T is the temperature, κ is the thermal conductivity, λ is the latent heat of evaporation/condensation, and Mw is the molecular weight of water. The effective volumetric heat capacity Cvt is defined by(3)Cvt=ϵCmgC+1-ϵCvf, where Cmg is the molar heat capacity and Cvf is the volumetric heat capacity of fiber.By Darcy’s law, the gas velocityug is defined as(4)ug=-kμ∂P∂x,where k is the permeability and μ is the dynamic viscosity, which usually is density-dependent for the compressible flow. Here we choose a linear form of μ≔νC, where ν is a certain constant.By the Hertz-Knudsen equation [11], the phase change rate Γce is defined as(5)Γce=βΓP-PsatTT,where βΓ is a positive constant, the saturation pressure Psat is determined from experimental measurements [12], and the pressure P is given by P=RCT, where R is the universal gas constant.With nondimensionalization, the sweat transport process (1)-(2) can be described by the following system:(6)ct-cθxx=-Γc,θ,0<x<L,t>0,(7)cθt+σθt-cθxθx-κθxx=λΓc,θ,0<x<L,t>0,where (·)x=∂/∂x, (·)t=∂/∂t, Γ(c,θ)=cθ-ps(θ), and ps(θ)~Psat(θ)/θ is a smooth and increasing function satisfying ps(0)=0.Since the right boundary is exposed to environment and the left boundary is connected to the body, we consider commonly used Robin type boundary conditions(8)cθx=α1c-μ1,x=0,0<t≤T,(9)cθx=α2μ2-c,x=L,0<t≤T,(10)θx=β1θ-ν1,x=0,0<t≤T,(11)θx=β2ν2-θ,x=L,0<t≤T,and the initial conditions(12)cx,0=c0x,θx,0=θ0x,0≤x≤L.Physically, parameters αi, βi, μi, νi, i=1,2, and σ are nonnegative constants [1, 2, 6]. We define initial condition parameters c0(x)≥c_,θ0(x)≥θ_ with c_ and θ_ being positive constants.Due to the strong nonlinearity and the coupling of the system, both theoretical and numerical analyses of the system are difficult. Numerical analysis for some related systems of parabolic/elliptic equations can be found in [13–20]. Existence and uniqueness of a classical solution for a steady-state model was given in [10]. Existence of a weak solution for the corresponding dynamic models was given in [21, 22]. Positivity of temperature and nonnegativity of vapor density were also proved here. Recently, a finite difference method second-order in space and first-order in time for the system (6)-(12) was presented in [23], where the backward semi-implicit Euler scheme is applied in the temporal direction and central finite difference approximations are used in the spatial direction. In [23], authors presented optimal error estimates under the assumption that the step size τ and h are smaller than a positive constant.In this paper, we propose an uncoupled leap-frog finite difference method for the system (6)-(12) with second-order accuracy in both spatial and temporal directions. We prove the existence and uniqueness of a solution to the finite difference system with optimal error estimates in the discrete L2 norm, under the condition that the mesh size τ and h are smaller than a positive constant which depends solely upon the physical parameters involved in the equations. Due to the strong nonlinearity and the coupling of equations, the method presented in [23] does not apply to the leap-frog scheme directly. One of the difficulties is to show convergence of the numerical solution without restriction on the grid ratio. In this paper, we assume that the solution (c(x,t),θ(x,t)) to the system (6)-(12) satisfies that (13)cmin≤cx,t≤cmax,θmin≤θx,t≤θmaxfor some positive constants cmin, cmax, θmax, and θmin.The manuscript is organized as follows: in Section2, we present an uncoupled leap-frog finite difference method for the nonlinear sweat transport system. In Section 3, we prove the existence and uniqueness of the solution to the sweat transport system with the optimal error estimate in the discrete L2 norm. Numerical results will be presented in Section 4 to support our theoretical results.


## 2. The Leap-Frog Finite Difference Scheme

For convenience of calculations, we add the equation (6) times θ into the equation (7); thus the governing system (6)-(7) can be modified as (14)ct-cθxx=-Γc,θ,0<x<L,t>0,(15)c+σθt-cθxθx-κθxx=λ+θΓc,θ,0<x<L,t>0. Due to the practical interest in a long time period, say 8–24 hours, we present an uncoupled leap-frog finite difference scheme in the temporal direction and the central finite difference (volume) scheme in the spatial direction for the above system with the initial/boundary conditions (8)-(12).LetT be a positive number, let Ωh={xi∣xi=ih,0≤i≤M} be a uniform partition in [0,L], and let Ωτ={tn∣tn=nτ,0≤n≤N} be a uniform partition in [0,T], where h=L/M and τ=T/N are the step size in the spatial and temporal directions, respectively. Denote xi+1/2=(xi+xi+1)/2 in the spatial cell and let {vin∣0≤i≤M,0≤n≤N} be a mesh function defined on Ωhτ, where Ωhτ=Ωh×Ωτ. Some notations are introduced below: (16)∇tvin=12τvin+1-vin-1,vin^=12vin+1+vin-1,vi+1/2n^=12vi+1n+1+vin-1,vi+1/2n=12vin+vi+1n,δxvi+1/2n=1hvi+1n-vin,δx∗vin^=1hvi+1/2n^-vi-1/2n^,δxvi+1/2n^=12δxvi+1/2n+1+δxvi+1/2n-1,δx2vi+1/2n^=1hδxvi+1/2n^-δxvi-1/2n^,from which(17)δxvi+1/2mvi+1/2n+vi+1/2mδxvi+1/2n=δxvi+1/2mvi+1/2n. The discrete system is defined by(18)∇tcin-1hci+1/2nδxθi+1/2n+θi+1/2nδxci+1/2n^-ci-1/2nδxθi-1/2n+θi-1/2nδxci-1/2n^=-Γcin,θin,1≤i≤M-1,(19)∇tc0n-2hc1/2nδxθ1/2n+θ1/2nδxc1/2n^-α1c0n^-μ1=-Γc0n,θ0n,(20)∇tcMn-2hα2μ2-cMn^-cM-1/2nδxθM-1/2n+θM-1/2nδxcM-1/2n^=-ΓcMn,θMn,(21)cin^+σ∇tθin-12ui-1/2n^δxθi-1/2n^+ui+1/2n^δxθi+1/2n^-κhδxθi+1/2n^-δxθi-1/2n^=λ+θinΓcin,θin,1≤i≤M-1,(22)c0n^+σ∇tθ0n-u1/2n^δxθ1/2n^-2κhδxθ1/2n^-β1θ0n^-ν1=λ+θ0nΓc0n,θ0n,(23)cMn^+σ∇tθMn-uM-1/2n^δxθM-1/2n^-2κhβ2ν2-θMn^-δxθM-1/2n^=λ+θMnΓcMn,θMn,1≤n≤N-1,and the discrete initial conditions(24)ci0=c0xi,θi0=θ0xi,(25)ci1=c0xi+τctxi,0,θi1=θ0xi+τθtxi,0,0≤i≤M,where(26)ui+1/2n^=ci+1/2nδxθi+1/2n+θi+1/2nδxci+1/2n^.The computational procedure of the uncoupled leap-frog scheme at each time step is listed below:Step 1.
The vapor concentrationcjn+1 can be calculated by solving the tridiagonal linear systems defined in (18)-(20).Step 2.
With the updated vapor concentrationcjn+1, we can get c0n^ and ui+1/2n^ correspondingly.Step 3.
Finally, the temperatureθjn+1 can be obtained by solving the tridiagonal linear system (21)-(23).


## 3. The Leap-Frog Scheme and the Optimal Error Estimate

In this section, we will show the existence and uniqueness of the solution to the system (18)-(26) with optimal error estimates in the discrete L2 norm. Let v={vjn}j=0M and z={zj}j=0M be two mesh functions on Ωh. We define the inner product and norms by (27)v,z=h12v0z0+∑i=1M-1vizi+12vMzM,v=h12v02+∑i=1M-1vi2+12vM2,v∞=max0≤i≤Mvi,δxv=h∑i=0M-1vi+1-vih2,δx2v=h∑i=1M-1δx2vi2.Let(C,Θ) be the solution of the system (6)-(12) and Cin=c(xi,tn), Θin=θ(xi,tn). The error functions are defined by(28)c~in=Cin-cin,θ~in=Θin-θin,0≤i≤M,0≤n≤N.We state our main result in the theorem below:Theorem 1.
Suppose that the solution(c,v) of the system (6)-(12) is in C4,3([0,L]×[0,T]), satisfying (13). Then there exist positive constants h0 and E0, independent of h and τ, such that, when τ≤Ech≤h0, the finite difference scheme (18)-(26) is uniquely solvable and(29)c~n+12+θ~n+12+τ∑m=1nδxc~m^2+δxθ~m^2≤E0τ2+h22,1≤n≤N-1.To prove the theorem, we make a stronger assumption that there existsE0>0, independent of n,h,τ, such that the inequality,(30)c~n+12+θ~n+12+δxθ~n+12≤E0τ2+h22,holds for n≤k-1. We prove the assumption and the theorem by induction method. By the initial condition (26), this is true for n=0. In the next subsection we will show that this is also true for n=1. In this part, we let Ec be a generic positive constant, which is associated with the physical parameters σ, λ, κ, cmin, cmax, θmin, θmax, the parameters involved in initial and boundary conditions and the solution of the system (6)-(12). Ec is independent of time step n, mesh size h, τ, and constant E0.


### 3.1. The Leap-Frog Scheme and Preliminaries

For convenience of calculations, we further introduce some notations. Letu=cθx,w=θx; thus the sweat transport system (6)-(7) can be reduced to(31)ct-ux=-Γc,θ,0≤x≤L,0<t≤T,(32)u=cθx+θcx,0≤x≤L,0<t≤T,(33)c+σθt-uθx-κwx=λ+θΓc,θ,0≤x≤L,0<t≤T,(34)w=θx,0≤x≤L,0<t≤T,with the initial and boundary conditions(35)u=α1c-μ1,x=0,0<t≤T,u=α2μ2-c,x=L,0<t≤T,w=β1θ-ν1,x=0,0<t≤T,w=β2ν2-θ,x=L,0<t≤T,cx,0=c0x,θx,0=θ0x,0≤x≤L.The discrete leap-frog system (18)-(23) is modified as(36)∇tcin-δx∗uin^=-Γcin,θin,1≤i≤M-1,∇tc0n-2hu1/2n^-α1c0n^-μ1=-Γc0n,θ0n,∇tcMn-2hα2μ2-cMn^-uM-1/2n^=-ΓcMn,θMn,ui+1/2n^=ci+1/2nδxθi+1/2n+θi+1/2nδxci+1/2n^,0≤i≤M-1,cin^+σ∇tθin-12ui-1/2n^δxθi-1/2n^+ui+1/2n^δxθi+1/2n^-κδx∗win^=λ+θinΓcin,θin,1≤i≤M-1,c0n^+σ∇tθ0n-u1/2n^δxθ1/2n^-2κhw1/2n^-β1θ0n^-ν1=λ+θ0nΓc0n,θ0n,cMn^+σ∇tθMn-uM-1/2n^δxθM-1/2n^-2κhβ2ν2-θMn^-wM-1/2n^=λ+θMnΓcMn,θMn,wi+1/2n^=δxθi+1/2n^,0≤i≤M-1,1≤n≤N-1.LetUin=u(xi,tn)=(cθ)x(xi,tn) and Win=w(xi,tn)=θx(xi,tn). We denote by uin and win the corresponding finite difference solution and(37)u~i+1/2n=Ui+1/2n-ui+1/2n,w~i+1/2n=Wi+1/2n-wi+1/2n,0≤i≤M-1,1≤n≤N-1.We get(38)∇tCin-δx∗Uin^=-ΓCin,Θin+Rc,in,1≤i≤M-1,∇tC0n-2hU1/2n^-α1C0n^-μ1=-ΓC0n,Θ0n+Rc,0n∇tCMn-2hα2μ2-CMn^-UM-1/2n^=-ΓCMn,ΘMn+Rc,Mn,Ui+1/2n^=Ci+1/2nδxΘi+1/2n+Θi+1/2nδxCi+1/2n^+Ru,i+1/2n,0≤i≤M-1,Cin^+σ∇tΘin-12Ui-1/2n^δxΘi-1/2n^+Ui+1/2n^δxΘi+1/2n^-κδx∗Win^=λ+ΘinΓCin,Θin+Rθ,in,1≤i≤M-1,C0n^+σ∇tΘ0n-U1/2nδxΘ1/2n^-2κhW1/2n^-β1Θ0n^-ν1=λ+Θ0nΓC0n,Θ0n+Rθ,0n,CMn^+σ∇tΘMn-UM-1/2n^δxΘM-1/2n^-2κhβ2ν2-ΘMn^-WM-1/2n^=λ+ΘMnΓCMn,ΘMn+Rθ,Mn,Wi+1/2n^=δxΘi+1/2n^+Rw,i+1/2n,0≤i≤M-1,1≤n≤N-1,and the initial conditions(39)Ci0=c0xi,Θi0=θ0xi,Ci1=c0xi+τctxi,0+τRc,i1,Θi1=θ0xi+τθtxi,0+τRθ,i1where (40)Rθ,i1,Rθ,i1≤E1τ+h2,0≤i≤M.Rc,in,Rθ,in≤E2τ2+h2,1≤i≤M-1,Ru,i+1/2n,Rw,i+1/2n≤E2τ2+h2,0≤i≤M-1,Rc,0n,Rc,Mn≤E2τ2+h2,Rθ,0n,Rθ,Mn≤E2τ2+h,1≤n≤N-1.Subtracting the system (36) from the system (38), we get the error equations(41)∇tc~in-δx∗u~in^=-ΓCin,Θin-Γcin,θin+Rc,in≔R¯c,in,1≤i≤M-1,(42)∇tc~0n-2hu~1/2n^-α1c~0n^=-ΓC0n,Θ0n-Γc0n,θ0n+Rc,0n≔R¯c,0n,(43)∇tc~Mn-2h-α2c~Mn^-u~M-1/2n^=-ΓCMn,ΘMn-ΓcMn,θMn+Rc,Mn≔R¯c,Mn,(44)u~i+1/2n^=c~i+1/2nδxΘi+1/2n+θ~i+1/2nδxCi+1/2n^+ci+1/2nδxθ~i+1/2n+θi+1/2nδxc~i+1/2n^+Ru,i+1/2n,0≤i≤M-1,(45)cin^+σ∇tθ~in-12ui-1/2n^δxθ~i-1/2n^+ui+1/2n^δxθ~i+1/2n^-κδx∗w~in^=-c~in^∇tΘin+12u~i-1/2n^δxΘi-1/2n^+u~i+1/2n^δxΘi+1/2n^+λ+ΘinΓCin,Θin-λ+θinΓcin,θin+Rθ,in≔R¯θ,in,1≤i≤M-1,(46)c0n^+σ∇tθ~0n-u1/2n^δxθ~1/2n^-2κhw~1/2n^-β1θ~0n^=-c~0n^∇tΘ0n^+u~1/2n^δxΘ1/2n^+λ+Θ0nΓC0n,Θ0n-λ+θ0nΓc0n,θ0n+Rθ,0n≔R¯θ,0n,(47)cMn^+σ∇tθ~Mn-uM-1/2n^δxθ~M-1/2n^-2κh-β2θ~Mn^-w~M-1/2n^=-c~Mn^∇tΘMn+u~M-1/2n^δxΘM-1/2n^+λ+ΘMnΓCMn,ΘMn-λ+θMnΓcMn,θMn+Rθ,Mn≔R¯θ,Mn,(48)w~i+1/2n^=δxθ~i+1/2n^+Rw,i+1/2n,0≤i≤M-1,1≤n≤N-1,and(49)c~i0=0,θ~i0=0,c~i1=τRc,i1,θ~i1=τRθ,i1,and by (40), we can directly derive the inequality(50)c~12+θ~12+δxθ~12≤E0τ2+h22. To prove our main theorem, the following formula will be often used:(51)v1/2z0+h∑i=1M-1δxvizi-vM-1/2zM=∑i=0M-1vi+1/2zik-zi+1k.In the following lemma, we present discrete Sobolev interpolation formulas and the proof can be found in [24].Lemma 2.
Letv and z be two mesh functions. Then for any positive constant ϵ,(52)v∞2≤ϵδxv2+1ϵ+1Lv2,(53)δxv∞≤ϵδx2v+Ecδxv2.Lemma 3.
(54) c ~ n ∞ , θ ~ n ∞ ≤ 3 E 0 1 / 2 τ 7 / 4 + h 3 / 2 ,  1 ≤ n ≤ k , (55) c ~ x n ^ , θ ~ x n ^ ≤ E 0 1 / 2 τ 3 / 2 + h ,  1 ≤ n ≤ k - 1 .Proof.
From (30) for 0≤n≤k-1, we have(56)c~n+12+τδxc~n+12≤E0τ2+h22,θ~n+12+τδxθ~n+12≤E0τ2+h22.When τ≤h, with the inverse inequality we have(57)c~in+12≤c~n+1∞2≤2h-1c~n+12≤2h-1E0τ2+h22≤8E0h3.When h≤τ, by taking ϵ=τ1/2 in Lemma 2,(58)c~in+12≤c~n+1∞2≤τ1/2δxc~n+12+τ-1/2+L-1c~n+12≤τ-1/2+τ1/2L-1E0τ2+h22≤8E0τ7/2,1≤i≤M-1,0≤n≤k-1.The first part of (54) is obtained and the second part and the inequality (55) can be proved similarly.In addition, by Lemma3, there exist constants E3>0 and s0>0 such that, when h,τ≤s0,(59)cmin2≤cin≤2cmax,θmin2≤θin≤2θmax,0≤i≤M-1,1≤n≤k,and(60)∇tΘin,∇tCin≤E3,0≤i≤M-1,1≤n≤k,δxCi+1/2n,δxΘi+1/2n,δxCi+1/2n^,δxΘi+1/2n^≤E3,0≤i≤M-1,1≤n≤k,ci+1/2n,θi+1/2n,≤E3,0≤i≤M-1,1≤n≤k,Γci+1/2n,θi+1/2n≤E3,0≤i≤M-1,1≤n≤k,λ+ΘinΓCin,Θin-λ+θinΓcin,θin≤E3c~in+θ~in,0≤i≤M,1≤n≤k,ΓCin,Θin-Γcin,θin≤E3c~in+θ~in,0≤i≤M,1≤n≤k.


### 3.2. The Existence and Uniqueness

Since the coefficient matrix in the system (18)-(20) is strictly diagonally dominant, thus the system (18)-(20) has a unique solution cik+1. Here we will discuss the boundedness of cik+1.Multiplying (41)-(43) by hc~in^, hc~0n^/2, and hc~Mn^/2, respectively, we get(61)∇tc~n,c~n^+α1c~0n^2+α2c~Mn^2=u~1/2n^c~0n^+h∑i=1M-1δx∗u~in^c~in^+u~M-1/2n^c~Mn^+h∑i=1M-1Rc,inc~in^+h2Rc,0nc~0n^+h2Rc,Mnc~Mn^-h12ΓC0n,Θ0n-Γc0n,θ0nc~0n^+∑i=1M-1ΓCin,Θin-Γcin,θinc~in^+12ΓCMn,ΘMn-ΓcMn,θMnc~Mn^≔J1+J2+J3,1≤n≤k.With (44), (51), (59), (60), and (40), we see that(62)-J1=-u~1/2n^c~0n^-h∑i=1M-1δxu~in^c~in^+u~M-1/2n^c~Mn^=h∑i=0M-1u~i+1/2n^δxc~i+1/2n^=h∑i=0M-1c~i+1/2nδxΘi+1/2n+θ~i+1/2nδxCi+1/2n^+ci+1/2nδxθ~i+1/2n+θi+1/2nδxc~i+1/2n^+Ru,i+1/2nδxc~i+1/2n^≥θmin2δxc~n^2-2cmaxδxθ~nδxc~n^-E3c~n+θ~nδxc~n^-h∑i=0M-1Ru,i+1/2nδxc~i+1/2n^≥θmin4δxc~n^2-Ecδxθ~n2+c~n2+θ~n2-Ecτ2+h22, and by using (60) again, we have(63)J3≤E3c~n+θ~nc~n^, and with (52),(64)J2≤12h∑i=1M-1Rc,in2+c~in^2+h22Rc,0n2+12c~0n^2+h22Rc,Mn2+12c~Mn^2≤c~n^∞2+Ecc~n^2+E2L2τ2+h22+E2h2τ2+h22≤θmin8δxc~n^2+Ecc~n^2+Ecτ2+h22. Substituting the last three equations into (61) results in(65)12∇tc~n2+θmin8δxc~n^2+α1c~0n^2+α2c~Mn^2≤Ecδxθ~n2+c~n2+θ~n2+c~n^2+Ecτ2+h22,1≤n≤k, where we have noted (∇tc~n,c~n^)=(1/4τ)(c~n+12-c~n-12)=(1/2)∇tc~n2. Moreover, by the assumption of the induction,(66)c~k+12+4τθmin8δxc~k^2+α1c~0k^2+α2c~Mk^2≤c~k-12+4Ecτδxθ~k2+c~k2+θ~k2+c~k^2+4τEcτ2+h22.Since we have the fact that c~k^2≤(1/2)(c~k+12+c~k-12), thus(67)1-2τEcc~k+12+τθmin2δxc~k^2≤E0+4τEc+4EcE0+6τEcE0τ2+h22. When τEc<1/4, we can get the inequality as(68)c~k+12+τθminδxc~k^2≤E4τ2+h22.Since E4 are independent of k, by (13) when h and τ are small enough,(69)cik+1≥0,0≤i≤M.Now we try to prove our main theorem. By noting (44), (60), (40), and Lemma 3,(70)u~n^2=∑i=0M-1hci+1/2nδxθ~i+1/2n+θi+1/2nδxc~i+1/2n^+c~i+1/2nδxΘi+1/2n+θ~i+1/2nδxCi+1/2n+Ru,i+1/2n2≤5∑i=0M-1h4cmax2δxθ~i+1/2n2+4θmax2δxc~i+1/2n^2+E32c~i+1/2n2+E32θ~i+1/2n2+Ru,i+1/2n2≤20cmax2δxθ~n2+20θmax2δxc~n^2+5E32c~n2+θ~n2+5E22Lτ2+h22,1≤n≤kWe can see that, when τ≤h, the assumption of induction and (68) show that(71)u~k^2≤40θmax2c~k+12h2+Ech4≤160θmax2E4h2+Ech4, and when h≤τ, by (68),(72)u~k^2≤20θmax2E4τ2+h22τθmin+Ecτ4≤80θmax2θminE4τ3+Ecτ4, which means there exists an E5 independent of k, such that(73)u~k^≤E5τ3/2+h.Multiplying the error equation (41) by δx∗uik^ leads to(74)δx∗u~k^2≤∑i=1M-1∇tc~ikδx∗uik^+∑i=1M-1R¯c,iδx∗uik^≤4∇tc~k2+8E32θ~k2+c~k2+34δx∗u~k^2+4E22Lτ2+h22,that is,(75)δx∗u~k^2≤16∇tc~k2+32E32θ~k2+c~k2+16E22Lτ2+h22.We can see that, when τ≤h,(76)u~k^∞≤h-1/2u~k^≤E5h+h1/2, and when h≤τ, with Lemma 2,(77)u~k^∞2≤δx∗u~k^2+1+L-1u~k^2≤E6τ2, where E6 is independent of k. Then there exists s0>0, when h,τ≤s0,(78)uk^∞=max1≤i≤Mui-1/2k^≤2E3.With a time step condition τ≤Ech, we can see that the coefficient matrix of the system (21)-(23) is strictly diagonally dominant. Thus this system has a unique solution θik+1.


### 3.3. The Optimal Error Estimate

We have proved the existence and uniqueness of the solution to the system and have derived the estimate (65) for c~n+1. In this part, we try to derive an estimate for θ~n+1.Multiplying (45)-(47) by hθ~in^, hθ~0n^/2, and hθ~Mn^/2, respectively, we try to estimate each term below:(79)h12c0n^+σ∇tθ~0nθ~0n^+∑i=1M-1cin^+σ∇tθ~inθ~in^+12cMn^+σ∇tθ~Mnθ~Mn^-h2u1/2n^δxθ~1/2n^θ~0n^+∑i=1M-1ui-1/2n^δxθ~i-1/2n^+ui+1/2n^δxθ~i+1/2n^θ~in^+uM-1/2n^δxθ~M-1/2n^θ~Mn^-κw~1/2n^θ~0n^+h∑i=1M-1δx∗w~i+1/2n^θ~in^-w~M-1/2n^θ~Mn^+κβ1θ~0n^2+κβ2θ~Mn^2=-h12c~0n^∇tΘ0nθ~0n^+∑i=1M-1c~in^∇tΘinθ~in^+12c~Mn^∇tΘMnθ~Mn^+h2u~1/2n^δxΘ1/2n^θ~0n^+∑i=1M-1u~i-1/2n^δxΘi-1/2n^+u~i+1/2n^δxΘi+1/2n^θ~in^+u~M-1/2n^δxΘM-1/2n^θ~Mn^+h2λ+Θ0nΓC0n,Θ0n-λ+θ0nΓc0n,θ0nθ~0n^+h∑i=1M-1λ+ΘinΓCin,Θin-λ+θinΓcin,θinθ~in^+h2λ+ΘMnΓCMn,ΘMn-λ+θMnΓcMn,θMnθ~Mn^+h∑i=1M-1Rθ,inθ~in^+h2Rθ,0nθ~0n^+h2Rθ,Mnθ~Mn^≔J3+J4+J5+J6,1≤n≤k.According to Lemma2, (40), (51), and (60), three terms on the left can be bounded by(80)h12c0n^+σ∇tθ~0nθ~0n^+∑i=1M-1cin^+σ∇tθ~inθ~in^+12cMn^+σ∇tθ~Mnθ~Mn^=h212c0n^+σ∇tθ~0n2+∑i=1M-1cin^+σ∇tθ~in2+12cMn^+σ∇tθ~Mn2,-h2u1/2n^δxθ~1/2n^θ~0n^+∑i=1M-1ui-1/2n^δxθ~i-1/2n^+ui+1/2n^δxθ~i+1/2n^θ~in^+uM-1/2n^δxθ~M-1/2n^θ~Mn^=-h∑i=0M-1ui+1/2n^δxθ~i+1/2n^θ~i+1/2n^≤2E3h∑i=0M-1δxθ~i+1/2n^θ~i+1/2n^≤κ4δxθ~n^2+Ecθ~n^2,and(81)-κw~1/2n^θ~0n^+h∑i=1M-1δx∗w~in^θ~in^-w~M-1/2n^θ~Mn^=κh∑i=0M-1δxθ~i+1/2n^2+κh∑i=0M-1Rw,i+1/2δxθ~i+1/2n^≥κδxθ~n^2-Ech∑i=0M-1Rw,i+1/2n2-κ4δxθ~n^2≥3κ4δxθ~n^2-Ecτ2+h22,1≤n≤k. By (70), for those terms in the right hand side, we obtain(82)J3≤E3c~n^2+θ~n^2,J4≤hE3u~1/2n^·θ~0n^+∑i=1M-1u~i-1/2n^+u~i+1/2n^θ~in^+u~M-1/2n^·θ~Mn^≤2E3θ~n^2+∑i=0M-1hu~i+1/2n^2≤40E3θmax2δxc~n^2+Ecδxθ~n2+θ~n+cn2+E2τ2+h22,J5≤E3c~n+θ~nθ~n^.J6≤12h∑i=1M-1Rθ,in2+θ~in^2+h22Rθ,0n2+12θ~0n^2+h22Rθ,Mn2+12θ~Mn^2≤θ~n^∞2+Ecθ~n^2+E2L2τ2+h22+E2h2τ2+h2≤κ4δxθ~n^2+θ~n^2+Ecτ2+h22.Taking the last six equations into (79), we obtain(83)h212c0n^+σ∇tθ~0n2+∑i=1M-1cin^+σ∇tθ~in2+12cMn^+σ∇tθ~Mn2+κ4δxθ~n^2+κβ1θ~0n^2+κβ2θ~Mn^2≤40E3θmax2δxc~n^2+Ecδxθ~n2+c~n2+θ~n2+c~n^2+θ~n^2+Ecτ2+h22,1≤n≤k.Besides, we introduce a notation asθ~i∗=(1/2)θ~in+12+θ~in-12, and adding the first three equations into (36) by hθ~i∗/2, hθ~0∗/4, and hθ~M∗/4, respectively, we have(84)h212∇tc0nθ~0∗+∑i=1M-1∇tcinθ~i∗+12∇tcMnθ~M∗-12u1/2n^θ~0∗+h∑i=1M-1δx∗ui+1/2n^θ~i∗-uM-1/2n^θ~M∗+α12c0n^-μ1θ~0∗+α22cMn^-μ2θ~M∗=-h212θ~0∗Γc0n,θ0n+∑i=1M-1θ~i∗Γcin,θin+12θ~M∗ΓcMn,θMn,1≤n≤k.We now estimate the terms in (84). By (51), we denote(85)J7≔-12u1/2n^θ~0∗+h∑i=1M-1δx∗ui+1/2n^θ~i∗-uM-1/2n^θ~M∗=12∑i=0M-1ui+1/2n^θ~i+1∗-θ~i∗≤E3h∑i=0M-1θ~i+1/2n+1δxθ~i+1/2n+1+θ~i+1/2n-1δxθ~i+1/2n-1≤E32δxθ~n+12+δxθ~n-12+θ~n+12+θ~n-12. Using (60) again, we get(86)-h212θ~0∗Γc0n,θ0n+∑i=1M-1θ~i∗Γcin,θin+12θ~M∗ΓcMn,θMn≤E34θ~n+12+θ~n-12, and with Lemma 2, we have(87)h212∇tc0nθ~0∗+∑i=1M-1∇tcinθ~i∗+12∇tcMnθ~M∗≤α1μ12θ~0∗+α2μ22θ~M∗+Ecθ~n+12+θ~n-12+δxθ~n+12+δxθ~n-12≤α1μ1+α2μ22θ~n+1∞2+θ~n-1∞2+Ecθ~n+12+θ~n-12+δxθ~n+12+δxθ~n-12≤Ecθ~n+12+θ~n-12+δxθ~n+12+δxθ~n-12.Moreover, by noting the fact that(88)cin^+σ∇tθ~in2+θ~i∗∇tcin=∇tcin+σθ~in2, adding (87) in (84) and using Lemma 2 again, we further get(89)h212∇tc0n+σθ~0n2+∑i=1M-1∇tcin+σθ~in2+12∇tcMn+σθ~Mn2+κ4δxθ~n^2+κβ1θ~0n^2+κβ2θ~Mn^2≤40E3θmax2δxc~n^2+Ecδxθ~n2+θ~n+12+c~n+12+θ~n2+c~n2+θ~n-12+c~n-12+Ecτ2+h22,1≤n≤k.Multiplying the last equation with s1=θmin/320E3θmax2 and adding it into equation (65), we get(90)∇tc~n2+θmin8δxc~n^2+2α1c~0n^2+2α2c~Mn^2+s1κβ1θ~0n^2+s1κβ2θ~Mn^2+s1h212∇tc0n+σθ~0n2+∑i=1M-1∇tcin+σθ~in2+12∇tcMn+σθ~Mn2+κs14δxθ~n^2≤Ecδxθ~n2+θ~n+12+c~n+12+θ~n2+c~n2+θ~n-12+c~n-12+Ecτ2+h22,1≤n≤k.Finally we estimateδxθ~n. Multiplying the error equation (45) by -hδx2θ~n^/(cin^+σ) and summing up the resulting equations for i=1,2,…,M-1, we have(91)-h∑i=1M-1∇tθ~inδx2θ~in^+κcin^+σh∑i=1M-1δx∗w~in^δx2θ~in^=hcin^+σ∑i=1M-1c~in^∇tΘinδx2θ~in^-h2cin^+σ∑i=1M-1u~i-1/2n^δxΘi-1/2n^+u~i+1/2n^δxΘi+1/2n^δx2θ~in^-hcin^+σ∑i=1M-1λ+ΘinΓCin,Θin-λ+θinΓcin,θin+Rθ,inδx2θ~in^-h2cin^+σ∑i=1M-1ui-1/2n^δxθ~i-1/2n^+ui+1/2n^δxθ~i+1/2n^δx2θ~in^,1≤n≤k.For the first term, we have(92)-h∑i=1M-1∇tθ~inδx2θ~in^=-h2τ∑i=1M-1θ~in+1δx2θ~in^+h2τ∑i=1M-1θ~in-1δx2θ~in^=-12τ∑i=1M-1θ~in+1δxθ~i+1/2n^-δxθ~i-1/2n^+12τ∑i=1M-1θ~in-1δxθ~i+1/2n^-δxθ~i-1/2n^.By (51), we have(93)-h∑i=1M-1∇tθ~inδx2θ~in^=-12τθ~Mn+1δxθ~M-1/2n^-θ~1n+1δxθ~1/2n^-∑i=1M-1δxθ~i+1/2n^θ~i+1n+1-θ~in+1+12τθ~Mn-1δxθ~M-1/2n^-θ~1n-1δxθ~1/2n^-∑i=1M-1δxθ~i+1/2n^θ~i+1n-1-θ~i+1n-1=-∇tθ~Mnδxθ~M-1/2n^+∇tθ~1nδxθ~1/2n^+14τh∑i=1M-1δxθ~i+1/2n+12-h∑i=1M-1δxθ~i+1/2n-12=12∇tδxθ~n2-∇tθ~Mnδxθ~M-1/2n^+∇tθ~0nδxθ~1/2n^. For the second term, we have(94)κcin^+σh∑i=1M-1δx∗w~in^δx2θ~in^≥κcmax+σh∑i=1M-1δx2θ~in^2+κcmax+σh∑i=1M-1δx∗Rw,inδx2θ~in^≥κ2cmax+σh∑i=1M-1δx2θ~in^2-Ecτ2+h22,where we noted the fact that δx∗Rw,in≤Ec(τ2+h2). From (94) we can get(95)12∇tδxθ~n2+κ2cmax+σδx2θ~n^2+∇tθ~0nδxθ~1/2n^-∇tθ~Mnδxθ~M-1/2n^≤hcin^+σ∑i=1M-1c~in^∇tΘinδx2θ~in^-h2cin^+σ∑i=1M-1u~i-1/2n^δxΘi-1/2n^+u~i+1/2n^δxΘi+1/2n^δx2θ~in^-hcin^+σ∑i=1M-1λ+ΘinΓCin,Θin-λ+θinΓcin,θin+Rθ,inδx2θ~in^-h2cin^+σ∑i=1M-1ui-1/2n^δxθ~i-1/2n^+ui+1/2n^δxθ~i+1/2n^δx2θ~in^+Ecτ2+h22≔J8+J9+J10+J11,1≤n≤k.Then we estimate the term ∇tθ~0nδxθ~1/2n^, -∇tθ~Mnδxθ~M-1/2n^ and J8 to J11, respectively. From (46), we have(96)δxθ~1/2n^=β1θ~0n^+h2κc0n^+σ∇tθ~0n-u1/2n^δxθ~1/2n^-R¯θ,0n-Rw,1/2n.A straightforward calculation with Lemma2 leads to(97)∇tθ~0nδxθ~1/2n^=β1∇tθ~0nθ~0n^+h2κc0n^+σ∇tθ~0n-u1/2n^δxθ~1/2n^-R¯θ,0n∇tθ~0n-Rw,1/2n∇tθ~0n≥β12∇tθ~0n2+hσ2κ∇tθ~0n2-h2κu1/2n^δxθ~1/2n^∇tθ~0n-h2κR¯θ,0n∇tθ~0n-∇tRw,1/2nθ~0n+Rw,1/2n+1-Rw,1/2n2τθ~0n+1+Rw,1/2n-Rw,1/2n-12τθ~0n-1≥β12∇tθ~0n2+hσ4κ∇tθ~0n2-h2κσu1/2n^δxθ~1/2n^2+R¯θ,0n2-∇tRw,1/2nθ~0n-Ecτ2+h2θ~n+1∞+θ~n-1∞≥β12∇tθ~0n2+hσ4κ∇tθ~0n2-∇tRw,1/2nθ~0n-2E32κσδxθ~n^2-50E32θmax2κσδxc~n^2-Ecθ~n+12+θ~n-12+δxθ~n+12+δxθ~n2+δxθ~n-12+c~n2+θ~n2+c~n^2-Ecτ2+h22,where we noted the fact that (Rw,1/2n+1-Rw,1/2n)/2τ≤Ec(τ2+h2) and (Rw,1/2n-Rw,1/2n-1)/2τ≤Ec(τ2+h2). Similarly we estimate -∇tθ~Mnδxθ~M-1/2n^. From (47), we have(98)-δxθ~M-1/2n^=β2θ~Mn^+h2κcMn^+σ∇tθ~0n-uM-1/2n^δxθ~M-1/2n^-R¯θ,Mn+Rw,M-1/2n, such that(99)-∇tθ~Mnδxθ~M-1/2n^≥β22∇tθ~Mn2+hσ4κ∇tθ~Mn2+∇tRw,M-1/2nθ~0n-2E32κσδxθ~n^2-50E32θmax2κσδxc~n^2-Ecτ2+h22-Ecθ~n+12+θ~n-12+δxθ~n+12+δxθ~n2+δxθ~n-12+c~n2+θ~n2+c~n^2.For those terms in the right hand side, with (70), we have(100)J8≤Ecc~n^2+κ16cmax+σδx2θ~n^2,(101)J9≤E32cmin+σh∑i=1M-1u~i-1/2n^δx2θ~in^+u~i+1/2n^δx2θ~in^≤κ16cmax+σδx2θ~n^2+16E32cmax+σcmin+σ2κh∑i=1M-1u~i+1/2n^2≤κ16cmax+σδx2θ~n^2+320E32cmax+σθmax2cmin+σ2κδxc~n^2+Ecδxθ~n2+c~n2+θ~n2+Ecτ2+h22,(102)J10≤κ16cmax+σδx2θ~n^2+Ecc~n2+θ~n2+Ecτ2+h22,(103)J11≤E3h∑i=1M-1δxθ~i-1/2n^δx2θ~i-1/2n^+δxθ~i+1/2n^δx2θ~i+1/2n^≤2E3h∑i=0M-132E3cmax+σκδxθ~i+1/2n^κ32E3cmax+σδx2θ~i+1/2n^≤κ16cmax+σδx2θ~n^2+64E32cmax+σκδxθ~n^2.Plugging the last six equations into (95), we get(104)12∇tδxθ~n2+κ4cmax+σδx2θ~n^2+β12∇tθ~0n2+β22∇tθ~Mn2+∇tRw,M-1/2nθ~Mn-Rw,1/2nθ~0n≤E7δxc~n^2+E8δxθ~n^2+Ecτ2+h22+Ecδxθ~n2+c~n2+c~n^2+θ~n2+θ~n+12+θ~n-12+δxθ~n+12+δxθ~n-12.Multiplying the last equation with 0≤s2≤E9 and adding it into (90), we have(105)s22∇tδxθ~n2+s2κ4cmax+σδx2θ~n^2+θmin16δxc~n^2+κs18δxθ~n^2+∇tc~n2+β1s22∇tθ~0n2+β2s22∇tθ~Mn2+s2∇t-Rw,1/2nθ~0n+Rw,M-1/2nθ~Mn+s1h212∇tc0n+σθ~0n2+∑i=1M-1∇tcin+σθ~in2+12∇tcMn+σθ~Mn2≤Ecδxθ~n+12+c~n+12+θ~n+12+Ecδxθ~n2+δxθ~n-12+c~n2+c~n-12+θ~n2+θ~n-12.+Ecτ2+h22.Lettings=min⁡(s2/2,s2κ/2(cmax+σ),θmin/8,κs1/4,1,s1,s2β1/2,s2β2/2), we get(106)s∇tδxθ~n2+s2δx2θ~n^2+s2δxc~n^2+s2δxθ~n^2+s∇tc~n2+s∇tθ~0n2+s∇tθ~Mn2+s∇t-Rw,1/2nθ~0n+Rw,M-1/2nθ~Mn+sh212∇tc0n+σθ~0n2+∑i=1M-1∇tcin+σθ~in2+12∇tcMn+σθ~Mn2≤Ecδxθ~n+12+δxθ~n2+δxθ~n-12+c~n+12+c~n2+c~n-12+θ~n+12+θ~n2+θ~n-12+Ecθ~0n+12+θ~0n2+θ~0n-12+θ~Mn+12+θ~Mn2+θ~Mn-12+Ecτ2+h22.Letting(107)Fn+1=sδxθ~n+12+sc~n+12+θ~0n2+θ~Mn2+sh212c0n+1+σθ~0n+12+∑i=1M-1cin+1+σθ~in+12+12cMn+1+σθ~Mn+12, we have(108)Fn+1-Fn-1+sτδxc~n^2+δxθ~n^2+δx2θ~n^2≤2τEcFn+1+Fn+Fn-1+2τEcτ2+h22, from which we can find α=(-τEc+1-3(τEc)2)/(1-2τEc) and λ=(τEc+1-3(τEc)2)/(1-2τEc), such that(109)Fn+1+αFn+sτ1-2τEcδxc~n^2+δxθ~n^2+δx2θ~n^2≤λFn+αFn-1+2τEc1-2τEcτ2+h22. When τEc<1/4,(110)Fn+1+αFn+∑m=1nτsδxc~m^2+δxθ~m^2+δx2θ~n^2≤e6EcTF1+F0+τ2+h22, and(111)sδxθ~n+12+sc~n+12+sh212c0n+1+σθ~0n+12+∑i=1M-1cin+1+σθ~in+12+12cMn+1+σθ~Mn+12+∑m=1nτsδxc~m^2+δxθ~m^2+δx2θ~n^2≤e6EcTτ2+h22. Since Ec is independent of E0, with(112)E0=e6EcTs, we find that (30) holds for k=n. The induction and the proof of the theorem are completed.


## 3.1. The Leap-Frog Scheme and Preliminaries

For convenience of calculations, we further introduce some notations. Letu=cθx,w=θx; thus the sweat transport system (6)-(7) can be reduced to(31)ct-ux=-Γc,θ,0≤x≤L,0<t≤T,(32)u=cθx+θcx,0≤x≤L,0<t≤T,(33)c+σθt-uθx-κwx=λ+θΓc,θ,0≤x≤L,0<t≤T,(34)w=θx,0≤x≤L,0<t≤T,with the initial and boundary conditions(35)u=α1c-μ1,x=0,0<t≤T,u=α2μ2-c,x=L,0<t≤T,w=β1θ-ν1,x=0,0<t≤T,w=β2ν2-θ,x=L,0<t≤T,cx,0=c0x,θx,0=θ0x,0≤x≤L.The discrete leap-frog system (18)-(23) is modified as(36)∇tcin-δx∗uin^=-Γcin,θin,1≤i≤M-1,∇tc0n-2hu1/2n^-α1c0n^-μ1=-Γc0n,θ0n,∇tcMn-2hα2μ2-cMn^-uM-1/2n^=-ΓcMn,θMn,ui+1/2n^=ci+1/2nδxθi+1/2n+θi+1/2nδxci+1/2n^,0≤i≤M-1,cin^+σ∇tθin-12ui-1/2n^δxθi-1/2n^+ui+1/2n^δxθi+1/2n^-κδx∗win^=λ+θinΓcin,θin,1≤i≤M-1,c0n^+σ∇tθ0n-u1/2n^δxθ1/2n^-2κhw1/2n^-β1θ0n^-ν1=λ+θ0nΓc0n,θ0n,cMn^+σ∇tθMn-uM-1/2n^δxθM-1/2n^-2κhβ2ν2-θMn^-wM-1/2n^=λ+θMnΓcMn,θMn,wi+1/2n^=δxθi+1/2n^,0≤i≤M-1,1≤n≤N-1.LetUin=u(xi,tn)=(cθ)x(xi,tn) and Win=w(xi,tn)=θx(xi,tn). We denote by uin and win the corresponding finite difference solution and(37)u~i+1/2n=Ui+1/2n-ui+1/2n,w~i+1/2n=Wi+1/2n-wi+1/2n,0≤i≤M-1,1≤n≤N-1.We get(38)∇tCin-δx∗Uin^=-ΓCin,Θin+Rc,in,1≤i≤M-1,∇tC0n-2hU1/2n^-α1C0n^-μ1=-ΓC0n,Θ0n+Rc,0n∇tCMn-2hα2μ2-CMn^-UM-1/2n^=-ΓCMn,ΘMn+Rc,Mn,Ui+1/2n^=Ci+1/2nδxΘi+1/2n+Θi+1/2nδxCi+1/2n^+Ru,i+1/2n,0≤i≤M-1,Cin^+σ∇tΘin-12Ui-1/2n^δxΘi-1/2n^+Ui+1/2n^δxΘi+1/2n^-κδx∗Win^=λ+ΘinΓCin,Θin+Rθ,in,1≤i≤M-1,C0n^+σ∇tΘ0n-U1/2nδxΘ1/2n^-2κhW1/2n^-β1Θ0n^-ν1=λ+Θ0nΓC0n,Θ0n+Rθ,0n,CMn^+σ∇tΘMn-UM-1/2n^δxΘM-1/2n^-2κhβ2ν2-ΘMn^-WM-1/2n^=λ+ΘMnΓCMn,ΘMn+Rθ,Mn,Wi+1/2n^=δxΘi+1/2n^+Rw,i+1/2n,0≤i≤M-1,1≤n≤N-1,and the initial conditions(39)Ci0=c0xi,Θi0=θ0xi,Ci1=c0xi+τctxi,0+τRc,i1,Θi1=θ0xi+τθtxi,0+τRθ,i1where (40)Rθ,i1,Rθ,i1≤E1τ+h2,0≤i≤M.Rc,in,Rθ,in≤E2τ2+h2,1≤i≤M-1,Ru,i+1/2n,Rw,i+1/2n≤E2τ2+h2,0≤i≤M-1,Rc,0n,Rc,Mn≤E2τ2+h2,Rθ,0n,Rθ,Mn≤E2τ2+h,1≤n≤N-1.Subtracting the system (36) from the system (38), we get the error equations(41)∇tc~in-δx∗u~in^=-ΓCin,Θin-Γcin,θin+Rc,in≔R¯c,in,1≤i≤M-1,(42)∇tc~0n-2hu~1/2n^-α1c~0n^=-ΓC0n,Θ0n-Γc0n,θ0n+Rc,0n≔R¯c,0n,(43)∇tc~Mn-2h-α2c~Mn^-u~M-1/2n^=-ΓCMn,ΘMn-ΓcMn,θMn+Rc,Mn≔R¯c,Mn,(44)u~i+1/2n^=c~i+1/2nδxΘi+1/2n+θ~i+1/2nδxCi+1/2n^+ci+1/2nδxθ~i+1/2n+θi+1/2nδxc~i+1/2n^+Ru,i+1/2n,0≤i≤M-1,(45)cin^+σ∇tθ~in-12ui-1/2n^δxθ~i-1/2n^+ui+1/2n^δxθ~i+1/2n^-κδx∗w~in^=-c~in^∇tΘin+12u~i-1/2n^δxΘi-1/2n^+u~i+1/2n^δxΘi+1/2n^+λ+ΘinΓCin,Θin-λ+θinΓcin,θin+Rθ,in≔R¯θ,in,1≤i≤M-1,(46)c0n^+σ∇tθ~0n-u1/2n^δxθ~1/2n^-2κhw~1/2n^-β1θ~0n^=-c~0n^∇tΘ0n^+u~1/2n^δxΘ1/2n^+λ+Θ0nΓC0n,Θ0n-λ+θ0nΓc0n,θ0n+Rθ,0n≔R¯θ,0n,(47)cMn^+σ∇tθ~Mn-uM-1/2n^δxθ~M-1/2n^-2κh-β2θ~Mn^-w~M-1/2n^=-c~Mn^∇tΘMn+u~M-1/2n^δxΘM-1/2n^+λ+ΘMnΓCMn,ΘMn-λ+θMnΓcMn,θMn+Rθ,Mn≔R¯θ,Mn,(48)w~i+1/2n^=δxθ~i+1/2n^+Rw,i+1/2n,0≤i≤M-1,1≤n≤N-1,and(49)c~i0=0,θ~i0=0,c~i1=τRc,i1,θ~i1=τRθ,i1,and by (40), we can directly derive the inequality(50)c~12+θ~12+δxθ~12≤E0τ2+h22. To prove our main theorem, the following formula will be often used:(51)v1/2z0+h∑i=1M-1δxvizi-vM-1/2zM=∑i=0M-1vi+1/2zik-zi+1k.In the following lemma, we present discrete Sobolev interpolation formulas and the proof can be found in [24].Lemma 2.
Letv and z be two mesh functions. Then for any positive constant ϵ,(52)v∞2≤ϵδxv2+1ϵ+1Lv2,(53)δxv∞≤ϵδx2v+Ecδxv2.Lemma 3.
(54) c ~ n ∞ , θ ~ n ∞ ≤ 3 E 0 1 / 2 τ 7 / 4 + h 3 / 2 ,  1 ≤ n ≤ k , (55) c ~ x n ^ , θ ~ x n ^ ≤ E 0 1 / 2 τ 3 / 2 + h ,  1 ≤ n ≤ k - 1 .Proof.
From (30) for 0≤n≤k-1, we have(56)c~n+12+τδxc~n+12≤E0τ2+h22,θ~n+12+τδxθ~n+12≤E0τ2+h22.When τ≤h, with the inverse inequality we have(57)c~in+12≤c~n+1∞2≤2h-1c~n+12≤2h-1E0τ2+h22≤8E0h3.When h≤τ, by taking ϵ=τ1/2 in Lemma 2,(58)c~in+12≤c~n+1∞2≤τ1/2δxc~n+12+τ-1/2+L-1c~n+12≤τ-1/2+τ1/2L-1E0τ2+h22≤8E0τ7/2,1≤i≤M-1,0≤n≤k-1.The first part of (54) is obtained and the second part and the inequality (55) can be proved similarly.In addition, by Lemma3, there exist constants E3>0 and s0>0 such that, when h,τ≤s0,(59)cmin2≤cin≤2cmax,θmin2≤θin≤2θmax,0≤i≤M-1,1≤n≤k,and(60)∇tΘin,∇tCin≤E3,0≤i≤M-1,1≤n≤k,δxCi+1/2n,δxΘi+1/2n,δxCi+1/2n^,δxΘi+1/2n^≤E3,0≤i≤M-1,1≤n≤k,ci+1/2n,θi+1/2n,≤E3,0≤i≤M-1,1≤n≤k,Γci+1/2n,θi+1/2n≤E3,0≤i≤M-1,1≤n≤k,λ+ΘinΓCin,Θin-λ+θinΓcin,θin≤E3c~in+θ~in,0≤i≤M,1≤n≤k,ΓCin,Θin-Γcin,θin≤E3c~in+θ~in,0≤i≤M,1≤n≤k.


## 3.2. The Existence and Uniqueness

Since the coefficient matrix in the system (18)-(20) is strictly diagonally dominant, thus the system (18)-(20) has a unique solution cik+1. Here we will discuss the boundedness of cik+1.Multiplying (41)-(43) by hc~in^, hc~0n^/2, and hc~Mn^/2, respectively, we get(61)∇tc~n,c~n^+α1c~0n^2+α2c~Mn^2=u~1/2n^c~0n^+h∑i=1M-1δx∗u~in^c~in^+u~M-1/2n^c~Mn^+h∑i=1M-1Rc,inc~in^+h2Rc,0nc~0n^+h2Rc,Mnc~Mn^-h12ΓC0n,Θ0n-Γc0n,θ0nc~0n^+∑i=1M-1ΓCin,Θin-Γcin,θinc~in^+12ΓCMn,ΘMn-ΓcMn,θMnc~Mn^≔J1+J2+J3,1≤n≤k.With (44), (51), (59), (60), and (40), we see that(62)-J1=-u~1/2n^c~0n^-h∑i=1M-1δxu~in^c~in^+u~M-1/2n^c~Mn^=h∑i=0M-1u~i+1/2n^δxc~i+1/2n^=h∑i=0M-1c~i+1/2nδxΘi+1/2n+θ~i+1/2nδxCi+1/2n^+ci+1/2nδxθ~i+1/2n+θi+1/2nδxc~i+1/2n^+Ru,i+1/2nδxc~i+1/2n^≥θmin2δxc~n^2-2cmaxδxθ~nδxc~n^-E3c~n+θ~nδxc~n^-h∑i=0M-1Ru,i+1/2nδxc~i+1/2n^≥θmin4δxc~n^2-Ecδxθ~n2+c~n2+θ~n2-Ecτ2+h22, and by using (60) again, we have(63)J3≤E3c~n+θ~nc~n^, and with (52),(64)J2≤12h∑i=1M-1Rc,in2+c~in^2+h22Rc,0n2+12c~0n^2+h22Rc,Mn2+12c~Mn^2≤c~n^∞2+Ecc~n^2+E2L2τ2+h22+E2h2τ2+h22≤θmin8δxc~n^2+Ecc~n^2+Ecτ2+h22. Substituting the last three equations into (61) results in(65)12∇tc~n2+θmin8δxc~n^2+α1c~0n^2+α2c~Mn^2≤Ecδxθ~n2+c~n2+θ~n2+c~n^2+Ecτ2+h22,1≤n≤k, where we have noted (∇tc~n,c~n^)=(1/4τ)(c~n+12-c~n-12)=(1/2)∇tc~n2. Moreover, by the assumption of the induction,(66)c~k+12+4τθmin8δxc~k^2+α1c~0k^2+α2c~Mk^2≤c~k-12+4Ecτδxθ~k2+c~k2+θ~k2+c~k^2+4τEcτ2+h22.Since we have the fact that c~k^2≤(1/2)(c~k+12+c~k-12), thus(67)1-2τEcc~k+12+τθmin2δxc~k^2≤E0+4τEc+4EcE0+6τEcE0τ2+h22. When τEc<1/4, we can get the inequality as(68)c~k+12+τθminδxc~k^2≤E4τ2+h22.Since E4 are independent of k, by (13) when h and τ are small enough,(69)cik+1≥0,0≤i≤M.Now we try to prove our main theorem. By noting (44), (60), (40), and Lemma 3,(70)u~n^2=∑i=0M-1hci+1/2nδxθ~i+1/2n+θi+1/2nδxc~i+1/2n^+c~i+1/2nδxΘi+1/2n+θ~i+1/2nδxCi+1/2n+Ru,i+1/2n2≤5∑i=0M-1h4cmax2δxθ~i+1/2n2+4θmax2δxc~i+1/2n^2+E32c~i+1/2n2+E32θ~i+1/2n2+Ru,i+1/2n2≤20cmax2δxθ~n2+20θmax2δxc~n^2+5E32c~n2+θ~n2+5E22Lτ2+h22,1≤n≤kWe can see that, when τ≤h, the assumption of induction and (68) show that(71)u~k^2≤40θmax2c~k+12h2+Ech4≤160θmax2E4h2+Ech4, and when h≤τ, by (68),(72)u~k^2≤20θmax2E4τ2+h22τθmin+Ecτ4≤80θmax2θminE4τ3+Ecτ4, which means there exists an E5 independent of k, such that(73)u~k^≤E5τ3/2+h.Multiplying the error equation (41) by δx∗uik^ leads to(74)δx∗u~k^2≤∑i=1M-1∇tc~ikδx∗uik^+∑i=1M-1R¯c,iδx∗uik^≤4∇tc~k2+8E32θ~k2+c~k2+34δx∗u~k^2+4E22Lτ2+h22,that is,(75)δx∗u~k^2≤16∇tc~k2+32E32θ~k2+c~k2+16E22Lτ2+h22.We can see that, when τ≤h,(76)u~k^∞≤h-1/2u~k^≤E5h+h1/2, and when h≤τ, with Lemma 2,(77)u~k^∞2≤δx∗u~k^2+1+L-1u~k^2≤E6τ2, where E6 is independent of k. Then there exists s0>0, when h,τ≤s0,(78)uk^∞=max1≤i≤Mui-1/2k^≤2E3.With a time step condition τ≤Ech, we can see that the coefficient matrix of the system (21)-(23) is strictly diagonally dominant. Thus this system has a unique solution θik+1.


## 3.3. The Optimal Error Estimate

We have proved the existence and uniqueness of the solution to the system and have derived the estimate (65) for c~n+1. In this part, we try to derive an estimate for θ~n+1.Multiplying (45)-(47) by hθ~in^, hθ~0n^/2, and hθ~Mn^/2, respectively, we try to estimate each term below:(79)h12c0n^+σ∇tθ~0nθ~0n^+∑i=1M-1cin^+σ∇tθ~inθ~in^+12cMn^+σ∇tθ~Mnθ~Mn^-h2u1/2n^δxθ~1/2n^θ~0n^+∑i=1M-1ui-1/2n^δxθ~i-1/2n^+ui+1/2n^δxθ~i+1/2n^θ~in^+uM-1/2n^δxθ~M-1/2n^θ~Mn^-κw~1/2n^θ~0n^+h∑i=1M-1δx∗w~i+1/2n^θ~in^-w~M-1/2n^θ~Mn^+κβ1θ~0n^2+κβ2θ~Mn^2=-h12c~0n^∇tΘ0nθ~0n^+∑i=1M-1c~in^∇tΘinθ~in^+12c~Mn^∇tΘMnθ~Mn^+h2u~1/2n^δxΘ1/2n^θ~0n^+∑i=1M-1u~i-1/2n^δxΘi-1/2n^+u~i+1/2n^δxΘi+1/2n^θ~in^+u~M-1/2n^δxΘM-1/2n^θ~Mn^+h2λ+Θ0nΓC0n,Θ0n-λ+θ0nΓc0n,θ0nθ~0n^+h∑i=1M-1λ+ΘinΓCin,Θin-λ+θinΓcin,θinθ~in^+h2λ+ΘMnΓCMn,ΘMn-λ+θMnΓcMn,θMnθ~Mn^+h∑i=1M-1Rθ,inθ~in^+h2Rθ,0nθ~0n^+h2Rθ,Mnθ~Mn^≔J3+J4+J5+J6,1≤n≤k.According to Lemma2, (40), (51), and (60), three terms on the left can be bounded by(80)h12c0n^+σ∇tθ~0nθ~0n^+∑i=1M-1cin^+σ∇tθ~inθ~in^+12cMn^+σ∇tθ~Mnθ~Mn^=h212c0n^+σ∇tθ~0n2+∑i=1M-1cin^+σ∇tθ~in2+12cMn^+σ∇tθ~Mn2,-h2u1/2n^δxθ~1/2n^θ~0n^+∑i=1M-1ui-1/2n^δxθ~i-1/2n^+ui+1/2n^δxθ~i+1/2n^θ~in^+uM-1/2n^δxθ~M-1/2n^θ~Mn^=-h∑i=0M-1ui+1/2n^δxθ~i+1/2n^θ~i+1/2n^≤2E3h∑i=0M-1δxθ~i+1/2n^θ~i+1/2n^≤κ4δxθ~n^2+Ecθ~n^2,and(81)-κw~1/2n^θ~0n^+h∑i=1M-1δx∗w~in^θ~in^-w~M-1/2n^θ~Mn^=κh∑i=0M-1δxθ~i+1/2n^2+κh∑i=0M-1Rw,i+1/2δxθ~i+1/2n^≥κδxθ~n^2-Ech∑i=0M-1Rw,i+1/2n2-κ4δxθ~n^2≥3κ4δxθ~n^2-Ecτ2+h22,1≤n≤k. By (70), for those terms in the right hand side, we obtain(82)J3≤E3c~n^2+θ~n^2,J4≤hE3u~1/2n^·θ~0n^+∑i=1M-1u~i-1/2n^+u~i+1/2n^θ~in^+u~M-1/2n^·θ~Mn^≤2E3θ~n^2+∑i=0M-1hu~i+1/2n^2≤40E3θmax2δxc~n^2+Ecδxθ~n2+θ~n+cn2+E2τ2+h22,J5≤E3c~n+θ~nθ~n^.J6≤12h∑i=1M-1Rθ,in2+θ~in^2+h22Rθ,0n2+12θ~0n^2+h22Rθ,Mn2+12θ~Mn^2≤θ~n^∞2+Ecθ~n^2+E2L2τ2+h22+E2h2τ2+h2≤κ4δxθ~n^2+θ~n^2+Ecτ2+h22.Taking the last six equations into (79), we obtain(83)h212c0n^+σ∇tθ~0n2+∑i=1M-1cin^+σ∇tθ~in2+12cMn^+σ∇tθ~Mn2+κ4δxθ~n^2+κβ1θ~0n^2+κβ2θ~Mn^2≤40E3θmax2δxc~n^2+Ecδxθ~n2+c~n2+θ~n2+c~n^2+θ~n^2+Ecτ2+h22,1≤n≤k.Besides, we introduce a notation asθ~i∗=(1/2)θ~in+12+θ~in-12, and adding the first three equations into (36) by hθ~i∗/2, hθ~0∗/4, and hθ~M∗/4, respectively, we have(84)h212∇tc0nθ~0∗+∑i=1M-1∇tcinθ~i∗+12∇tcMnθ~M∗-12u1/2n^θ~0∗+h∑i=1M-1δx∗ui+1/2n^θ~i∗-uM-1/2n^θ~M∗+α12c0n^-μ1θ~0∗+α22cMn^-μ2θ~M∗=-h212θ~0∗Γc0n,θ0n+∑i=1M-1θ~i∗Γcin,θin+12θ~M∗ΓcMn,θMn,1≤n≤k.We now estimate the terms in (84). By (51), we denote(85)J7≔-12u1/2n^θ~0∗+h∑i=1M-1δx∗ui+1/2n^θ~i∗-uM-1/2n^θ~M∗=12∑i=0M-1ui+1/2n^θ~i+1∗-θ~i∗≤E3h∑i=0M-1θ~i+1/2n+1δxθ~i+1/2n+1+θ~i+1/2n-1δxθ~i+1/2n-1≤E32δxθ~n+12+δxθ~n-12+θ~n+12+θ~n-12. Using (60) again, we get(86)-h212θ~0∗Γc0n,θ0n+∑i=1M-1θ~i∗Γcin,θin+12θ~M∗ΓcMn,θMn≤E34θ~n+12+θ~n-12, and with Lemma 2, we have(87)h212∇tc0nθ~0∗+∑i=1M-1∇tcinθ~i∗+12∇tcMnθ~M∗≤α1μ12θ~0∗+α2μ22θ~M∗+Ecθ~n+12+θ~n-12+δxθ~n+12+δxθ~n-12≤α1μ1+α2μ22θ~n+1∞2+θ~n-1∞2+Ecθ~n+12+θ~n-12+δxθ~n+12+δxθ~n-12≤Ecθ~n+12+θ~n-12+δxθ~n+12+δxθ~n-12.Moreover, by noting the fact that(88)cin^+σ∇tθ~in2+θ~i∗∇tcin=∇tcin+σθ~in2, adding (87) in (84) and using Lemma 2 again, we further get(89)h212∇tc0n+σθ~0n2+∑i=1M-1∇tcin+σθ~in2+12∇tcMn+σθ~Mn2+κ4δxθ~n^2+κβ1θ~0n^2+κβ2θ~Mn^2≤40E3θmax2δxc~n^2+Ecδxθ~n2+θ~n+12+c~n+12+θ~n2+c~n2+θ~n-12+c~n-12+Ecτ2+h22,1≤n≤k.Multiplying the last equation with s1=θmin/320E3θmax2 and adding it into equation (65), we get(90)∇tc~n2+θmin8δxc~n^2+2α1c~0n^2+2α2c~Mn^2+s1κβ1θ~0n^2+s1κβ2θ~Mn^2+s1h212∇tc0n+σθ~0n2+∑i=1M-1∇tcin+σθ~in2+12∇tcMn+σθ~Mn2+κs14δxθ~n^2≤Ecδxθ~n2+θ~n+12+c~n+12+θ~n2+c~n2+θ~n-12+c~n-12+Ecτ2+h22,1≤n≤k.Finally we estimateδxθ~n. Multiplying the error equation (45) by -hδx2θ~n^/(cin^+σ) and summing up the resulting equations for i=1,2,…,M-1, we have(91)-h∑i=1M-1∇tθ~inδx2θ~in^+κcin^+σh∑i=1M-1δx∗w~in^δx2θ~in^=hcin^+σ∑i=1M-1c~in^∇tΘinδx2θ~in^-h2cin^+σ∑i=1M-1u~i-1/2n^δxΘi-1/2n^+u~i+1/2n^δxΘi+1/2n^δx2θ~in^-hcin^+σ∑i=1M-1λ+ΘinΓCin,Θin-λ+θinΓcin,θin+Rθ,inδx2θ~in^-h2cin^+σ∑i=1M-1ui-1/2n^δxθ~i-1/2n^+ui+1/2n^δxθ~i+1/2n^δx2θ~in^,1≤n≤k.For the first term, we have(92)-h∑i=1M-1∇tθ~inδx2θ~in^=-h2τ∑i=1M-1θ~in+1δx2θ~in^+h2τ∑i=1M-1θ~in-1δx2θ~in^=-12τ∑i=1M-1θ~in+1δxθ~i+1/2n^-δxθ~i-1/2n^+12τ∑i=1M-1θ~in-1δxθ~i+1/2n^-δxθ~i-1/2n^.By (51), we have(93)-h∑i=1M-1∇tθ~inδx2θ~in^=-12τθ~Mn+1δxθ~M-1/2n^-θ~1n+1δxθ~1/2n^-∑i=1M-1δxθ~i+1/2n^θ~i+1n+1-θ~in+1+12τθ~Mn-1δxθ~M-1/2n^-θ~1n-1δxθ~1/2n^-∑i=1M-1δxθ~i+1/2n^θ~i+1n-1-θ~i+1n-1=-∇tθ~Mnδxθ~M-1/2n^+∇tθ~1nδxθ~1/2n^+14τh∑i=1M-1δxθ~i+1/2n+12-h∑i=1M-1δxθ~i+1/2n-12=12∇tδxθ~n2-∇tθ~Mnδxθ~M-1/2n^+∇tθ~0nδxθ~1/2n^. For the second term, we have(94)κcin^+σh∑i=1M-1δx∗w~in^δx2θ~in^≥κcmax+σh∑i=1M-1δx2θ~in^2+κcmax+σh∑i=1M-1δx∗Rw,inδx2θ~in^≥κ2cmax+σh∑i=1M-1δx2θ~in^2-Ecτ2+h22,where we noted the fact that δx∗Rw,in≤Ec(τ2+h2). From (94) we can get(95)12∇tδxθ~n2+κ2cmax+σδx2θ~n^2+∇tθ~0nδxθ~1/2n^-∇tθ~Mnδxθ~M-1/2n^≤hcin^+σ∑i=1M-1c~in^∇tΘinδx2θ~in^-h2cin^+σ∑i=1M-1u~i-1/2n^δxΘi-1/2n^+u~i+1/2n^δxΘi+1/2n^δx2θ~in^-hcin^+σ∑i=1M-1λ+ΘinΓCin,Θin-λ+θinΓcin,θin+Rθ,inδx2θ~in^-h2cin^+σ∑i=1M-1ui-1/2n^δxθ~i-1/2n^+ui+1/2n^δxθ~i+1/2n^δx2θ~in^+Ecτ2+h22≔J8+J9+J10+J11,1≤n≤k.Then we estimate the term ∇tθ~0nδxθ~1/2n^, -∇tθ~Mnδxθ~M-1/2n^ and J8 to J11, respectively. From (46), we have(96)δxθ~1/2n^=β1θ~0n^+h2κc0n^+σ∇tθ~0n-u1/2n^δxθ~1/2n^-R¯θ,0n-Rw,1/2n.A straightforward calculation with Lemma2 leads to(97)∇tθ~0nδxθ~1/2n^=β1∇tθ~0nθ~0n^+h2κc0n^+σ∇tθ~0n-u1/2n^δxθ~1/2n^-R¯θ,0n∇tθ~0n-Rw,1/2n∇tθ~0n≥β12∇tθ~0n2+hσ2κ∇tθ~0n2-h2κu1/2n^δxθ~1/2n^∇tθ~0n-h2κR¯θ,0n∇tθ~0n-∇tRw,1/2nθ~0n+Rw,1/2n+1-Rw,1/2n2τθ~0n+1+Rw,1/2n-Rw,1/2n-12τθ~0n-1≥β12∇tθ~0n2+hσ4κ∇tθ~0n2-h2κσu1/2n^δxθ~1/2n^2+R¯θ,0n2-∇tRw,1/2nθ~0n-Ecτ2+h2θ~n+1∞+θ~n-1∞≥β12∇tθ~0n2+hσ4κ∇tθ~0n2-∇tRw,1/2nθ~0n-2E32κσδxθ~n^2-50E32θmax2κσδxc~n^2-Ecθ~n+12+θ~n-12+δxθ~n+12+δxθ~n2+δxθ~n-12+c~n2+θ~n2+c~n^2-Ecτ2+h22,where we noted the fact that (Rw,1/2n+1-Rw,1/2n)/2τ≤Ec(τ2+h2) and (Rw,1/2n-Rw,1/2n-1)/2τ≤Ec(τ2+h2). Similarly we estimate -∇tθ~Mnδxθ~M-1/2n^. From (47), we have(98)-δxθ~M-1/2n^=β2θ~Mn^+h2κcMn^+σ∇tθ~0n-uM-1/2n^δxθ~M-1/2n^-R¯θ,Mn+Rw,M-1/2n, such that(99)-∇tθ~Mnδxθ~M-1/2n^≥β22∇tθ~Mn2+hσ4κ∇tθ~Mn2+∇tRw,M-1/2nθ~0n-2E32κσδxθ~n^2-50E32θmax2κσδxc~n^2-Ecτ2+h22-Ecθ~n+12+θ~n-12+δxθ~n+12+δxθ~n2+δxθ~n-12+c~n2+θ~n2+c~n^2.For those terms in the right hand side, with (70), we have(100)J8≤Ecc~n^2+κ16cmax+σδx2θ~n^2,(101)J9≤E32cmin+σh∑i=1M-1u~i-1/2n^δx2θ~in^+u~i+1/2n^δx2θ~in^≤κ16cmax+σδx2θ~n^2+16E32cmax+σcmin+σ2κh∑i=1M-1u~i+1/2n^2≤κ16cmax+σδx2θ~n^2+320E32cmax+σθmax2cmin+σ2κδxc~n^2+Ecδxθ~n2+c~n2+θ~n2+Ecτ2+h22,(102)J10≤κ16cmax+σδx2θ~n^2+Ecc~n2+θ~n2+Ecτ2+h22,(103)J11≤E3h∑i=1M-1δxθ~i-1/2n^δx2θ~i-1/2n^+δxθ~i+1/2n^δx2θ~i+1/2n^≤2E3h∑i=0M-132E3cmax+σκδxθ~i+1/2n^κ32E3cmax+σδx2θ~i+1/2n^≤κ16cmax+σδx2θ~n^2+64E32cmax+σκδxθ~n^2.Plugging the last six equations into (95), we get(104)12∇tδxθ~n2+κ4cmax+σδx2θ~n^2+β12∇tθ~0n2+β22∇tθ~Mn2+∇tRw,M-1/2nθ~Mn-Rw,1/2nθ~0n≤E7δxc~n^2+E8δxθ~n^2+Ecτ2+h22+Ecδxθ~n2+c~n2+c~n^2+θ~n2+θ~n+12+θ~n-12+δxθ~n+12+δxθ~n-12.Multiplying the last equation with 0≤s2≤E9 and adding it into (90), we have(105)s22∇tδxθ~n2+s2κ4cmax+σδx2θ~n^2+θmin16δxc~n^2+κs18δxθ~n^2+∇tc~n2+β1s22∇tθ~0n2+β2s22∇tθ~Mn2+s2∇t-Rw,1/2nθ~0n+Rw,M-1/2nθ~Mn+s1h212∇tc0n+σθ~0n2+∑i=1M-1∇tcin+σθ~in2+12∇tcMn+σθ~Mn2≤Ecδxθ~n+12+c~n+12+θ~n+12+Ecδxθ~n2+δxθ~n-12+c~n2+c~n-12+θ~n2+θ~n-12.+Ecτ2+h22.Lettings=min⁡(s2/2,s2κ/2(cmax+σ),θmin/8,κs1/4,1,s1,s2β1/2,s2β2/2), we get(106)s∇tδxθ~n2+s2δx2θ~n^2+s2δxc~n^2+s2δxθ~n^2+s∇tc~n2+s∇tθ~0n2+s∇tθ~Mn2+s∇t-Rw,1/2nθ~0n+Rw,M-1/2nθ~Mn+sh212∇tc0n+σθ~0n2+∑i=1M-1∇tcin+σθ~in2+12∇tcMn+σθ~Mn2≤Ecδxθ~n+12+δxθ~n2+δxθ~n-12+c~n+12+c~n2+c~n-12+θ~n+12+θ~n2+θ~n-12+Ecθ~0n+12+θ~0n2+θ~0n-12+θ~Mn+12+θ~Mn2+θ~Mn-12+Ecτ2+h22.Letting(107)Fn+1=sδxθ~n+12+sc~n+12+θ~0n2+θ~Mn2+sh212c0n+1+σθ~0n+12+∑i=1M-1cin+1+σθ~in+12+12cMn+1+σθ~Mn+12, we have(108)Fn+1-Fn-1+sτδxc~n^2+δxθ~n^2+δx2θ~n^2≤2τEcFn+1+Fn+Fn-1+2τEcτ2+h22, from which we can find α=(-τEc+1-3(τEc)2)/(1-2τEc) and λ=(τEc+1-3(τEc)2)/(1-2τEc), such that(109)Fn+1+αFn+sτ1-2τEcδxc~n^2+δxθ~n^2+δx2θ~n^2≤λFn+αFn-1+2τEc1-2τEcτ2+h22. When τEc<1/4,(110)Fn+1+αFn+∑m=1nτsδxc~m^2+δxθ~m^2+δx2θ~n^2≤e6EcTF1+F0+τ2+h22, and(111)sδxθ~n+12+sc~n+12+sh212c0n+1+σθ~0n+12+∑i=1M-1cin+1+σθ~in+12+12cMn+1+σθ~Mn+12+∑m=1nτsδxc~m^2+δxθ~m^2+δx2θ~n^2≤e6EcTτ2+h22. Since Ec is independent of E0, with(112)E0=e6EcTs, we find that (30) holds for k=n. The induction and the proof of the theorem are completed.


## 4. Numerical Examples

We now numerically evaluate the performance of the proposed leap-frog scheme.Example 1.
First, we test the accuracy of our algorithm in an artificial example which is taken from [23]. The system is (113)Ct+ugCx=-Γ+fCx,t,C+σTt-κTxx+ugCxTx=λ+σTΓ+fTx,t,with the boundary conditions (8)-(12), where fC, fT, μi, and νi are coefficients decided by the exact solution(114)Cx,t=e-0.72tx21-x2+a1-x+bx,Tx,t=50e-0.72tx21-x2+c1-x+dx,with a,b,c,d being constants.
We apply the uncoupled leap-frog finite difference method to solve the artificial example. We chooseT=1 and L=1. Since the proposed scheme is of the second order in both spatial and temporal directions, we take τ=h such that the error bound is proportional to h2. We present the L2-norm errors and the order of convergence hr in Table 1 with h=L/200,L/400,L/800 at different time level. We can see clearly from Table 1 that the L2-norm errors for both components are proportional to h2, which confirms our theoretical analysis.Table 1
Numerical results of Example1 with τ=h.
c ~ n ^ t h = L / 200 h = L / 400 h=L/800 Orderr t = 0.25 2.914E-03 1.004E-03 0.254E-03 1.99 t = 0.50 8.163E-03 2.178E-03 0.527E-03 2.03 t = 0.75 1.397E-02 0.350E-02 0.084E-02 2.04 t = 1.00 2.032E-02 0.496E-02 0.118E-02 2.05  θ~n^ t h = L / 200 h = L / 400 h=L/800 Orderr  t = 0.25 9.509E-03 3.795E-03 0.947E-03 2.00 t = 0.50 3.860E-03 1.042E-03 0.253E-03 2.03 t = 0.75 7.115E-02 1.799E-02 0.439E-02 2.02 t = 1.00 1.069E-01 0.265E-01 0.065E-01 2.02Example 2.
In the second example, we discuss a typical clothing assembly in the textile industry [2, 4, 25]. The clothing assembly consists of three layers, in the middle is porous fibrous media, and the outside cover is exposed to a cold environment with fixed temperature and relative humidity while the inside cover is exposed to a mixture of air and vapor at higher temperature and relative humidity. In this paper, polyester porous media with laminated or nylon cover materials are tested. To compare with the experimental data in [12], a water equation is added to equations (1)-(2):(115)∂∂tρw1-ϵ′W=MwΓcewhere W is water content, ρw is the density of water, ϵ is the porosity with liquid water content, and ϵ′ is the porosity without liquid water content. We have(116)ϵ=ϵ′-ρfρwW1-ϵ′and the effective heat conductivity is defined by(117)κ=ϵκg+1-ϵκswhere κg is the thermal conductivity of gas and κs is the thermal conductivity of the fiber-water mixture [2, 6, 7], given by(118)κs=1ρf+Wρw1ρfκf+Wρwκw-1. The values of these physical parameters for polyester media are presented in Table 2. Other parameters values can be found in [2, 6, 7].
The initial conditions for the vapor, temperature, and water content are given by(119)C=65%PsatTRT,T=25∘C,W=0,att=0.
We apply the uncoupled leap-frog finite difference method for solving the sweat transport system defined in (1)-(2) coupled with the water equation (115). Since only the right side of the water equation includes c and θ, therefore, the water equation is calculated separately. Numerically, at each time step, we first find solution cjn+1, θjn+1 by procedure (18)-(26), and then Wjn+1 can be solved by following nonnormalized discrete formate:(120)ρf1-ϵ′jnτWjn+1-Wjn=MwΓce,jn. Then we evaluate the parameters explicitly in (18)-(26) based on Wjn+1. Here all numerical results are obtained by taking the time step size τ=20s and spatial mesh size h=L/100. We present numerical results of vapor, temperature, and water content at 8 hours and 24 hours, respectively, for the porous polyester media assembly with laminated cover in Figure 1 and with nylon cover in Figure 2. The comparisons between numerical results of water content and experimental measurements [12] are given in last two subfigures, where the blue lines represent the numerical solution and the red line is given by experimental measurement.Table 2
Physical parameters for batting materials.
parameter polyester unit ρf 1.39 × 1 0 3 kgm−3 ρw 1 × 1 0 3 kgm−3 ϵ′ 0.993 κf 1 × 1 0 - 1 W m−1K−1 κw 5.7 × 1 0 - 1 W m−1K−1 Cvf 1.17 × 1 0 6 J m−1K−1 L 4.92 × 1 0 - 2 mFigure 1
Numerical results for 10-pile polyester batting sandwiched by two layers of laminated fabric.Figure 2
Numerical results for 10-pile polyester batting sandwiched by two layers of nylon fabric.


## 5. Conclusion

As a subsequent work of [23], we have presented an uncoupled leap-frog finite difference method for the sweat transport system in porous textile media, which is governed by a strongly coupled, nonlinear parabolic system. Optimal L2 error estimates were presented, which imply that the numerical scheme is unconditionally stable. Both theoretical analysis and numerical example indicate that the current scheme is second order accurate in both the temporal and spatial directions. Since the scheme is decoupled for the system, the method can be applied efficiently for problems in higher-dimensional space. Under certain time-step restrictions, the analysis can also be extended to the multidimensional problems.



---

*Source: 8649308-2019-06-12.xml*