Research article Special Issues

Optimized Schwarz waveform relaxation for heterogeneous Cattaneo-Vernotte non-Fourier heat transfer

  • The non-Fourier heat transfer in heterogeneous media is crucial for material science and biomedical engineering. The optimized Schwarz waveform relaxation (OSWR) method is an efficient approach for solving such problems due to its divide-and-conquer strategy. Despite the wave-type nature of non-Fourier heat transfer, the short phase-lag time leads to more parabolic-like behavior. To address this, in the OSWR method, we employed Robin boundary conditions to transmit information along the interface. Using Fourier analysis, we derived and rigorously optimized the convergence factors of the OSWR algorithm with scaled Robin and Robin-Robin transmission conditions. The resulting optimized transmission parameters were provided in explicit form for direct application in the OSWR algorithm, along with corresponding convergence factor estimates. Interestingly, the results show that a larger heterogeneity contrast actually accelerates the convergence, rather than deteriorating it. Furthermore, the OSWR algorithm with the Robin-Robin condition exhibits mesh-independent convergence asymptotically. However, the presence of the phase-lag time is found to slow down the convergence of the OSWR algorithm. These theoretical findings were validated through numerical experiments.

    Citation: Feng Hu, Yingxiang Xu. Optimized Schwarz waveform relaxation for heterogeneous Cattaneo-Vernotte non-Fourier heat transfer[J]. AIMS Mathematics, 2025, 10(3): 7370-7395. doi: 10.3934/math.2025338

    Related Papers:

    [1] Afraz Hussain Majeed, Sadia Irshad, Bagh Ali, Ahmed Kadhim Hussein, Nehad Ali Shah, Thongchai Botmart . Numerical investigations of nonlinear Maxwell fluid flow in the presence of non-Fourier heat flux theory: Keller box-based simulations. AIMS Mathematics, 2023, 8(5): 12559-12575. doi: 10.3934/math.2023631
    [2] Abdelbaki Choucha, Asma Alharbi, Bahri Cherif, Rashid Jan, Salah Boulaaras . Decay rate of the solutions to the Bresse-Cattaneo system with distributed delay. AIMS Mathematics, 2023, 8(8): 17890-17913. doi: 10.3934/math.2023911
    [3] W. Abbas, Ahmed M. Megahed, Osama M. Morsy, M. A. Ibrahim, Ahmed A. M. Said . Dissipative Williamson fluid flow with double diffusive Cattaneo-Christov model due to a slippery stretching sheet embedded in a porous medium. AIMS Mathematics, 2022, 7(12): 20781-20796. doi: 10.3934/math.20221139
    [4] Junjiang Lai, Zhencheng Fan . Stability for discrete time waveform relaxation methods based on Euler schemes. AIMS Mathematics, 2023, 8(10): 23713-23733. doi: 10.3934/math.20231206
    [5] S. R. Mishra, Subhajit Panda, Mansoor Alshehri, Nehad Ali Shah, Jae Dong Chung . Sensitivity analysis on optimizing heat transfer rate in hybrid nanofluid flow over a permeable surface for the power law heat flux model: Response surface methodology with ANOVA test. AIMS Mathematics, 2024, 9(5): 12700-12725. doi: 10.3934/math.2024621
    [6] Danhua Wang, Kewang Chen . Decay properties for the Cauchy problem of the linear JMGT-viscoelastic plate with Cattaneo heat conduction. AIMS Mathematics, 2025, 10(5): 12079-12091. doi: 10.3934/math.2025547
    [7] Yudhveer Singh, Devendra Kumar, Kanak Modi, Vinod Gill . A new approach to solve Cattaneo-Hristov diffusion model and fractional diffusion equations with Hilfer-Prabhakar derivative. AIMS Mathematics, 2020, 5(2): 843-855. doi: 10.3934/math.2020057
    [8] Yousef Jawarneh, Humaira Yasmin, M. Mossa Al-Sawalha, Rasool Shah, Asfandyar Khan . Numerical analysis of fractional heat transfer and porous media equations within Caputo-Fabrizio operator. AIMS Mathematics, 2023, 8(11): 26543-26560. doi: 10.3934/math.20231356
    [9] Khalil Ur Rehman, Wasfi Shatanawi, Zeeshan Asghar, Haitham M. S. Bahaidarah . Neural networking analysis for MHD mixed convection Casson flow past a multiple surfaces: A numerical solution. AIMS Mathematics, 2023, 8(7): 15805-15823. doi: 10.3934/math.2023807
    [10] Khalid Masood . Exploring assisting and opposing flow characteristics in a saturated darcy medium with non-uniform heat source or sink and suspended microbes. AIMS Mathematics, 2025, 10(6): 14804-14839. doi: 10.3934/math.2025666
  • The non-Fourier heat transfer in heterogeneous media is crucial for material science and biomedical engineering. The optimized Schwarz waveform relaxation (OSWR) method is an efficient approach for solving such problems due to its divide-and-conquer strategy. Despite the wave-type nature of non-Fourier heat transfer, the short phase-lag time leads to more parabolic-like behavior. To address this, in the OSWR method, we employed Robin boundary conditions to transmit information along the interface. Using Fourier analysis, we derived and rigorously optimized the convergence factors of the OSWR algorithm with scaled Robin and Robin-Robin transmission conditions. The resulting optimized transmission parameters were provided in explicit form for direct application in the OSWR algorithm, along with corresponding convergence factor estimates. Interestingly, the results show that a larger heterogeneity contrast actually accelerates the convergence, rather than deteriorating it. Furthermore, the OSWR algorithm with the Robin-Robin condition exhibits mesh-independent convergence asymptotically. However, the presence of the phase-lag time is found to slow down the convergence of the OSWR algorithm. These theoretical findings were validated through numerical experiments.



    Non-Fourier heat transfer refers to an intricate heat transfer process that occurs under extreme conditions like high temperature or high pressure. This phenomenon deviates from the traditional Fourier heat conduction theory and plays important roles in diverse scientific fields including biomedicine [1,2] and materials science [3,4,5]. The classical Fourier heat law, expressed as

    q(x,t)=κTx,

    states that any heat flux generated at a specific location is instantaneously felt throughout the medium, where κ is the thermal conductivity. However, this behavior is inadequate in certain scenarios, such as heat transfer in micro-structures [3], because heat transfer involves irregular collisions between microscopic particles, atoms, or phonons, and the average communication time between them is non-negligible on this extraordinarily small time scale. Therefore, the Fourier law is insufficient to fully capture the real physics of this process. To solve the problem, Cattaneo [6] and Vernotte [7] independently proposed a modified heat flux principle,

    q(x,t+τ)=κTx, (1.1)

    which extends the Fourier law by introducing the phase-lag time τ in the relation of temperature gradient and heat flux. Physically, τ represents the finite buildup time for the commencement of heat flow in a medium. This implies that there is a time delay in the response of the heat flux to the temperature gradient, and the thermal propagation speed κ/τ is no longer infinite [8]. The phase-lag time τ can be determined by measuring the amplitude and phase as a function of the modulation frequency [9]. As pointed out in the literature, the phase-lag time τ is generally on a very small scale, especially in the application of nanocomposites, for example in microseconds or picoseconds [10,11]. Taylor expanding (1.1) in τ to the first order, and applying the energy conservation equation

    q(x,t)x=ρcTt,

    one arrives at the Cattaneo-Vernotte (CV) equation

    τ2Tt2+Tt=β2Tx2, (1.2)

    where β=κρc, and ρ and c denote the density and the heat capacity, respectively.

    In this paper, we are interested in the non-Fourier heat transfer in multi-layer media, where τ and β are discontinuous across the interface formed by adjacent media. To simplify the analysis, we consider only the case where the defining domain Ω is a bounded open domain in Rd (d=1,2,3) consisting of two layers occupied in Ω1 and Ω2, i.e., ¯Ω=¯Ω1¯Ω2 and Ω1Ω2=, and τ and β are piecewise constants in domain Ω,

    (τ,β)={(τ1,β1)inΩ1,(τ2,β2)inΩ2.

    The solution and the heat flux are continuous across the interface Γ=Ω1Ω2:

    [T]Γ=T1|ΓT2|Γ=0,[βTn]Γ=β1T1n1|Γ+β2T2n2|Γ=0,

    where Ti denotes the solution in subdomain Ωi, and ni is the normal vector pointing to outside of Ωi. To summarize, the CV model problem under consideration reads for i=1,2 and j=3i:

    τi2Tit2+Titβi2Tix2=fiinΩi×(0,T],Ti(x,0)=Ti0inΩi,Tit(x,0)=˜Ti0inΩi,Ti(x,t)Tj(x,t)=0onΓ×(0,T],βiTini+βjTjnj=0onΓ×(0,T]. (1.3)

    The above settings emerge in many applications, including biomedical treatment [12,13,14,15] and heat conduction of composite materials [16,17,18].

    It is crucial to numerically solve the heterogeneous CV Eq (1.3), as analytical solving is infeasible due to the complex coupling condition on the interface and/or the intricate defining domain. Two decades ago, Duhamel proposed a finite integral transform pair to solve the hyperbolic heat conduction problem in heterogeneous media [19]. Recently, solving heterogeneous problems similar to (1.3) has attracted significant attention. Singh et al. considered the second-order wave equation with piecewise discontinuous coefficients and proposed to solve it using a compact difference scheme, achieving second-order convergence in time and fourth-order convergence in space [20]. A hybrid numerical method based on the Haar wavelet and finite difference was used to solve the telegraph interface model with discontinuous coefficients [21]. Aiming to improve the performance near the interface, several works have focused on the (fitted) finite element method for 1D or 2D non-Fourier bioheat transfer models in multi-layered media, where optimal error estimates for finite element solutions in L2 and H1 norms have been obtained [22,23,24,25]. Additionally, the finite volume method has been investigated for solving the non-Fourier nonlinear heat conduction in a heterogeneous medium, where the thermal conductivity and the heat capacity are temperature-dependent [26].

    The CV model (1.3) is essentially a typical interface problem. Therefore, classical methods for interface problems, such as the fitted finite element method and the immersed boundary element method, can be applied, coupled with a promising time-stepping strategy. These methods form a large-scale algebraic system over the whole defining domain and generally encounter the problem of heavy ill-conditioning due to the parameter contrast between different media, which should be further tackled by preconditioning [27,28]. In addition, discretizing the space domain monolithically often leads to non-physical oscillations in the solution near the interface. Domain decomposition methods are adaptive to tackle this issue since only homogeneous subproblems should be solved iteratively. Among them, the optimized Schwarz method [29] attracted tremendous attention due to its excellent performance, since Robin-type transmission conditions are involved and optimized toward fast convergence. Maday and Magoulès extended the optimized Schwarz method to solve diffusion models in heterogeneous media, where a min-max problem was solved numerically to determine the optimal transmission parameters [30]. Later, Gander and Dubois [31] derived the explicit expression of the optimized transmission parameters in several different transmission conditions, as well as the corresponding convergence rate estimates, which helped to find that the strong heterogeneity between subdomains accelerates the convergence of the optimized Schwartz method. Similar results were also recovered by Gander and Vanzan for a more general problem: a diffusion-reaction diffusion coupling model [32]. Note that when the model problem evolves in time, the optimized Schwarz waveform relaxation (OSWR) method should be applied, which decomposes the large space-time problem into homogeneous subproblems and solves them iteratively via transmission conditions, even adapts to solving nonlinear problems [33]. For convection-diffusion problems with discontinuous coefficients, the OSWR method was analyzed in detail in [34], see also [35,36] for an OSWR coupled with discontinuous Galerkin time stepping, where an energy analysis was used to determine the optimized transmission parameters. We remark here that the OSWR algorithm for the homogeneous CV model has been rigorously analyzed under the name of the telegrapher equation [37].

    Motivated by the advantages of OSWR for solving heterogeneous problems, in this paper, we aim at developing an OSWR algorithm for solving the CV model (1.3), which leads to heat transfer subproblems on homogeneous media to be solved iteratively via well-designed transmission conditions. We focus only on the case where the phase-lag time τ1,2 is on a small scale, which corresponds to the non-Fourier heat transfer in composites consisting of multi-layered solid materials. The transmission parameters involved in the OSWR algorithm with scaled Robin and Robin-Robin transmission conditions are rigorously optimized. We note here that the Robin-Robin condition provides more degree of freedom to improve the algorithm's performance through solving a much harder optimization problem, while keeping its implementation and computational cost at each iteration the same as the scaled Robin condition. The corresponding convergence factor estimate reveals several good properties of the OSWR algorithm: the stronger heterogeneity leads to faster convergence for both transmission conditions, and the OSWR algorithm achieves mesh-independent convergence if the Robin-Robin condition is applied. We comment here that several parallel-in-time strategies such as the ParaDIAG [38], the parareal [39], or the unmapped tent pitching [40] algorithms can be applied to solve this subproblem and lead to more efficient computation; however, this is beyond the scope of this paper.

    The rest of the paper is organized as follows. In Section 2 we describe the OSWR algorithm for the CV model (1.3) and perform a Fourier analysis to find the convergence factor in the Fourier frequency domain. In Section 3, we propose using the scaled Robin condition to transfer information between subdomains and analyze the optimized transmission parameters. In Section 4, we consider a Robin-Robin transmission condition and analyze the optimized parameters, as well as the corresponding convergence factor estimate. The theoretical results are illustrated using numerical experiments in Section 5, and we conclude in Section 6.

    The defining domain Ω of the CV Eq (1.3) can be naturally decomposed into two subdomains Ω1 and Ω2. Specially, in each subdomain, only a homogeneous CV equation should be solved. Thus, we propose to solve the interface coupling problem using the OSWR method: for n=1,2,..., iteratively solve:

    τ12Tn1t2+Tn1tβ12Tn1x2=f1inΩ1×(0,T],Tn1(x,0)=T10inΩ1,Tn1t(x,0)=˜T10inΩ1,B1(Tn1,Tn12)=0 onΓ×(0,T], (2.1)

    and

    τ22Tn2t2+Tn2tβ22Tn2x2=f2inΩ2×(0,T],Tn2(x,0)=T20inΩ2,Tn2t(x,0)=˜T20inΩ2,B2(Tn2,Tn11)=0onΓ×(0,T], (2.2)

    where T0i,i=1,2, are initial values on the interface Γ×(0,T] with Γ=¯Ω1¯Ω2, and Bi,i=1,2, are transmission operators along the interface, which should be determined toward fast convergence of the algorithm. Noting that, since the phase-lag τi is on a very small scale, Eq (1.3) behaves more like a parabolic equation. Since the Robin conditions are effective for parabolic problems in transmitting information along the interface [41,42], we specify the transmission conditions B1(Tn1,Tn12) and B2(Tn2,Tn11) as the following Robin-type conditions:

    β1n1Tn1+S1Tn1=β2n2Tn12+S1Tn12onΓ×(0,T], (2.3)
    β2n2Tn2+S2Tn2=β1n1Tn11+S2Tn11onΓ×(0,T], (2.4)

    where Si,i=1,2, are operators on the interface. Once converged, we have in matrix form

    [IS1IS2][β1n1T1T1]=[IS1IS2][β2n2T2T2], (2.5)

    where the limit limnTni=Ti is used. Reformulating (2.5) as

    [IS1IS2]([β1n1T1T1][I00I][β2n2T2T2])=0

    shows that if S1 and S2 are chosen such that the coefficient matrix [IS1IS2] is invertible, the transmission conditions (2.3) and (2.4), once converged, will recover the interface condition in (1.3). Extending the solution by zero for t<0, the Fourier analysis* can be used to describe the convergence of the OSWR algorithms (2.1)–(2.4), and the operators S1 and S2 can then be determined toward fast convergence through optimization. Though the techniques also apply to high-dimensional problems, to well illustrate the analysis, we focus here only on a one-dimensional example. When extended to higher dimensions, a Fourier transform in multiple directions (e.g., two directions for two-dimensional problems) should be performed, resulting in a more complex optimization problem to be analyzed. We then assume that the subdomains are infinite one-dimensional domains: Ω1=(,0) and Ω2=(0,). Highlighted by [44], the infinite domain assumption only slightly affects the optimized transmission operators Si. Performing Fourier transforms with frequency k in the time direction on the governing equations in (2.1) and (2.2), respectively, gives, for kK,

    *Since the time t must be positive, it is natural to perform a Laplace transform with a complex frequency. However, justified by [43], these two transforms are equivalent in analyzing the OSWR method. Thus, for ease of analysis, the Fourier transform is adopted here.

    τ1k2ˆTn1(x,k)+ikˆTn1(x,k)β1xxˆTn1(x,k)=0x(,0),τ2k2ˆTn2(x,k)+ikˆTn2(x,k)β2xxˆTn2(x,k)=0x(0,), (2.6)

    where the homogeneous governing equation f=0 is considered, which means that we are analyzing the error equation, and K represents the set of Fourier frequencies involved in the calculation: K=KK+ with K=[kmax,kmin] and K+=[kmin,kmax], where kmin and kmax should be estimated according to the mesh size and the boundary conditions applied. Similarly, Fourier transforming the transmission conditions (2.3) and (2.4) gives, for kK,

    (β1x+σ1(k))ˆTn1(0,k)=(β2x+σ1(k))ˆTn12(0,k),(β2xσ2(k))ˆTn2(0,k)=(β1xσ2(k))ˆTn11(0,k), (2.7)

    where σi,i=1,2, are Fourier symbols of Si. Solving (2.6) one finds

    ˆTn1(x,k)=An1eλ1x,ˆTn2(x,k)=An2eλ2x, (2.8)

    where λi=τik2+ikβi and An1 and An2 depend on k. Inserting Eq (2.8) into Eq (2.7), we get

    An1=β2λ2+σ1β1λ1+σ1An12, An2=β1λ1σ2β2λ2σ2An11,

    which allow us to establish the following relation:

    Ani=ρ(k,β1,β2,τ1,τ2,σ1,σ2)An1i, i=1,2,

    where the convergence factor

    ρ(k,β1,β2,τ1,τ2,σ1,σ2)=β2λ2σ1β1λ1+σ1β1λ1σ2β2λ2+σ2.

    For ease of mathematical analysis, we would like to consider only the case τ:=τ1=τ2, which would greatly simplify the analysis while keeping the analytical result still applicable. We thus define

    λ:=β1λ1=β2λ2=τk2+ik,

    and the convergence factor simplifies to

    ρ(k,β1,β2,σ1,σ2)=β2λσ1β1λ+σ1β1λσ2β2λ+σ2.

    Choosing

    σ1=σopt1:=β2λ,σ2=σopt2:=β1λ (2.9)

    leads the convergence factor to ρ0 for all k, which means the OSWR algorithms (2.1) and (2.2) converge in two iterations. However, the corresponding transmission operators Si,i=1,2, are nonlocal and expensive to use, so we aim to find their local approximations that are easy to apply. We note here that if σi,i=1,2, are real numbers, it is easy to check that |ρ(k,β1,β2,σ1,σ2)| is an even function in k, thus one can approximate σi by a real constant pi and determine the optimal transmission parameters p1,p2 by minimizing the modulus of the convergence factor over all positive frequencies

    minp1,p2RmaxkK+|ρ(k,β1,β2,p1,p2)|.

    Guided by (2.9), we choose

    σ1(k)=β2p,σ2(k)=β1p, (3.1)

    where p is a constant to be determined by

    minpRmaxkK+ρR(k,p,v), (3.2)

    with

    ρR(k,p,v)=(ap)2+b2(a+pv)2+b2(a+pv)2+b2,

    v=β2β1 describes the coefficient contrast in the heterogeneous media, a and b are the real and imaginary parts of λ(k),

    a=2k2τk2+kτ2k2+1,b=τk2+kτ2k2+12=k2a.

    In the setting (3.1) the Robin constant p is scaled by β3i, hence the name "scaled Robin".

    Theorem 3.1. For positive p>0, ρR(k,p,v)<1 for all frequencies kK+, that is to say, the OSWR algorithms (2.1) and (2.2) converge for all positive p. In addition, the min-max problem (3.2) is equivalent to

    minp>0maxkK+ρR(k,p,v). (3.3)

    Proof. The proof of the first assertion is shown by the inequality: for p>0,

    0<ρR(k,p,v)<(ap)2+b2(a+p)2+b2<1.

    For the second assertion, it is easy to check that ρR(k,p,v)<ρR(k,p,v) for positive p>0. Moreover, we have ρR(k,0,v)1>ρR(k,p,v) for p>0. Thus, it is sufficient to consider the min-max problem (3.2) for positive p, i.e., the second assertion holds.

    Lemma 3.1 (Restricting the range of p). If p[a2min+b2min,a2max+b2max], then p is not a solution of the min-max problem (3.3), where amin=a(kmin), amax=a(kmax), bmin=b(kmin), and bmax=b(kmax).

    Proof. Highlighted by Theorem 3.1, we only need to consider the positive p>0. Taking the derivative of ρR in p, it shows

    ρRp=(va3+p(v2+1)a2+v(p2+b2)a+b2p(v1)2)(v+1)2(a2+b2p2)(p2v2+2apv+a2+b2)32((a2+b2)v2+2apv+p2)(a2+b2)v2+2apv+p2v2. (3.4)

    Since a,b, and v are all positive, we can deduce that ρRp<0 if p<a2min+b2min. Hence, ρR is a decreasing function in p, and cannot attain its minimum for p<a2min+b2min because decreasing p would increase ρR for all kK+. For p>a2max+b2max, a similar analysis can be performed, too. Thus we conclude that the solution of (3.3) must lie in the interval [a2min+b2min,a2max+b2max].

    Next, we need to analyze the maxima of ρR(k,p,v). We claim that ρR(k,p,v) attains its local maximum at either kmin or kmax. However, the complicated expression of ρR makes it impossible to analyze the maximum directly. We thus aim to analyze with the aid of an approximation of ρR. Using Taylor expansion we get, for small τ,

    a=2k224k32τ+O(τ2),b=2k2+24k32τ+O(τ2). (3.5)

    Then, replacing a and b in ρR with aapp=bapp=2k2 gives an approximate convergence factor

    ρappR(k,p,v)=(aappp)2+b2app(aapp+pv)2+b2app(aapp+pv)2+b2app.

    The following theorem shows the approximation error.

    Theorem 3.2. (Approximation error) There exists a constant C depending on the parameters k,p, and v but independent of τ such that the following estimate holds:

    |ρR(k,p,v)ρappR(k,p,v)|Cτ. (3.6)

    Proof. Using Taylor expansion in τ, one finds

    (ap)2+b2(a+pv)2+b2=(aappp)2+b2app(aapp+pv)2+b2app+2p(v+1)(p2v+k)k324p22kp+k(p2v2+2kpv+k)32τ+O(τ2), (3.7)

    and

    (ap)2+b2(a+pv)2+b2=(aappp)2+b2app(aapp+pv)2+b2app+2p(v+1)(p2+kv)k324p22kp+k(p2+2kpv+kv2)32vτ+O(τ2). (3.8)

    Multiplying (3.7) by (3.8), a simple calculation leads to the desired result (3.6).

    Lemma 3.2. The convergence factor ρR(k,p,v) achieves its maximum at either kmin or kmax in an asymptotic sense of τ small, more precisely,

    maxkK+ρR(k,p,v)=max{ρappR(kmin,p,v),ρappR(kmax,p,v)}+O(τ). (3.9)

    Proof. Taking the derivative of ρappR in k directly gives

    ρappRk=(v+1)2(2p(v2v+1)k322p3(v2v+1)k2(p4vk2v))p2kpv+kv2+p2v2k(p2v2+2kpv+k)32(4kv2+42kpv+4p2)).

    Investigating the signs of kρappR, we conclude that ρappR(k,p,v) decreases in k monotonically for k[kmin,p2] and increases in k monotonically for k[p2,kmax]. Hence ρappR must achieve its maximum at either kmin or kmax. Consequently, (3.9) follows from an application of Theorem 3.2.

    Recall that in the CV model (1.3), the phase-lag time τ is very small, and in the following analysis, we will regard kmin and kmax as the local maximum points of ρR.

    Theorem 3.3. (Solve by equi-oscillation) For small τ, the min-max problem (3.3) is solved by equi-oscillating the convergence factor ρR at the two endpoints kmin and kmax,

    ρR(kmin,p,v)=ρR(kmax,p,v). (3.10)

    In addition, the solution p is unique.

    Proof. Using Lemma 3.1 one finds that (3.3) is equivalent to the following min-max problem:

    minp[a2min+b2min,a2max+b2max](maxkK+ρR(k,p,v)). (3.11)

    For all p[a2min+b2min,a2max+b2max], using (3.4) we get

    ρRp(kmin,p,v)>0,ρRp(kmax,p,v)<0.

    Therefore, increasing p will increase ρR(kmin,p,v) and decrease ρR(kmax,p,v) for p[a2min+b2min,a2max+b2max]. Noting as well the following facts:

    ρR(kmin,a2min+b2min,v)<ρR(kmax,a2min+b2min,v),ρR(kmin,a2max+b2max,v)>ρR(kmax,a2max+b2max,v),

    we conclude that there exists a unique p satisfying the equi-oscillation Eq (3.10), and we thus arrive at our assertion.

    Theorem 3.4. (Optimized scaled Robin transmission conditions) For sufficiently large kmax, the unique solution p to the equi-oscillation problem (3.10) has the following asymptotic expression:

    p=214a12mink14max, (3.12)

    which leads the convergence factor ρR to the following estimate:

    maxkK+ρR(k,p,v)=1amin(v+1)2214vk14max+O(k12max). (3.13)

    Proof. Guided by numerical tests, we make the ansatz p=Cpk14max. Using asymptotic expansion for kmax large, one finds

    ρR(kmin,p,v)=1amin(v+1)2Cpvk14max+O(k12max) (3.14)

    and

    ρR(kmax,p,v)=12Cp(v+1)22vk14max+O(k12max).

    Matching the coefficients of k14max in the above two equations leads to the solution p as shown in (3.12). Inserting (3.12) into (3.14) we arrive at the desired result.

    Remark 3.1. The convergence factor estimate (3.13) reveals that a large coefficient contrast v leads to a larger coefficient of the term k14max, which ultimately results in a smaller convergence factor estimate. In other words, a stronger heterogeneity contrast leads to faster convergence of the OSWR algorithm. In addition, by taking a limit as v1, the result degenerates to that for CV equations defined on homogeneous media, i.e., β1=β2. This is a new result that has not been published publicly.

    Remark 3.2. The phase-lag τ is implicitly contained in a(kmin), and thus it enters both the optimized parameter (3.12) and the convergence factor estimate (3.13). Noting that the function a(k) decreases in τ, we conclude from (3.13) that a relatively large phase-lag τ will slow down the convergence of the OSWR algorithm. Specifically, for τ=0, the analysis degenerates to the case of a heterogeneous Fourier heat conduction problem, which is investigated in [45] using an equi-oscillation principle and numerical optimization. In contrast, we advanced the investigation using the asymptotic analysis technique, which allows us to establish the equivalence between the min-max problem and the equi-oscillation problem, tackling the challenge introduced by the phase-lag and providing the optimized parameter in explicit form that is ready to use in applications. More interestingly, the obtained convergence rate estimate allows us to reveal the effects of heterogeneity contrast, as shown in Remark 3.1.

    We now consider the following parameter setting:

    σ1(k)=β2p,σ2(k)=β1q,

    where p,q are constants to be determined by

    minp,qRmaxkK+ρRR(k,p,q,v), (4.1)

    with the convergence factor

    ρRR(k,p,q,v)=(ap)2+b2(a+pv)2+b2(aq)2+b2(a+qv)2+b2.

    In this case, the subdomain problems use different Robin transmission parameters, hence the name "Robin-Robin". Unlike solving by two-point equi-oscillation for the scaled Robin condition, since two free parameters are involved, the min-max problem (4.1) will be solved by three-point equi-oscillation, which is of course more challenging. To tackle this issue, we decompose the two-variable optimization problem into two sequential optimizations of single variable, and analyze them using asymptotic analysis together with an appropriate convergence factor approximation.

    Theorem 4.1. For all p,q>0, the convergence factor ρRR(k,p,q,v)<1 for all frequencies kK+, i.e., the OSWR algorithms (2.1) and (2.2) converge. In addition, the min-max problem (4.1) is equivalent to

    minp,q>0maxkK+ρRR(k,p,q,v). (4.2)

    Proof. The first assertion is justified by the following inequality:

    0<ρRR(k,p,q,v)<(ap)2+b2(a+p)2+b2(aq)2+b2(a+q)2+b2<1.

    Now we consider the second assertion. It is easy to check that for any fixed q and p>0, ρRR(k,p,q,v)<ρRR(k,p,q,v). While for any fixed p and q>0, ρRR(k,p,q,v)<ρRR(k,p,q,v). If both p and q are negative, we clearly have ρRR(k,p,q,v)>1. Moreover, because of ρRRp(k,0,q)<0 and ρRRq(k,p,0)<0, we find that (p,q)=(0,0) cannot solve (4.1). To sum up, the second assertion holds.

    Lemma 4.1. (Restricting the range of p and q) The solution (p,q) of the min-max problem (4.2) must be contained in [φ1(kmin),φ1(kmax)]×[φ2(kmin),φ2(kmax)], where φ1(k) and φ2(k) are defined as

    φ1(k)=(a2+b2)(v1)+(a2+b2)(a2(v+1)2+b2(v1)2)2av, (4.3)
    φ2(k)=(a2+b2)(1v)+(a2+b2)(a2(v+1)2+b2(v1)2)2a. (4.4)

    Proof. A direct calculation shows

    ρRRp=(aq)2+b2(v+1)(p(b2+a2ap)v+(a2+b2)(ap))((a+pv)2+b2)32(ap)2+b2(av+q)2+b2v2v2

    and the root of h(p):=p(a2+b2ap)v+(a2+b2)(ap) is φ1(k) defined in (4.3). Noting that for v>0, φ1(k)>0 for any k>kmin, thus for p<φ1(kmin), decreasing p increases the convergence factor ρRR for all k[kmin,kmax], and for p>φ1(kmax), increasing p increases the convergence factor. To conclude, the best p for minimizing ρRR must be contained in [φ1(kmin),φ1(kmax)]. A similar discussion shows that the best parameter q for minimizing ρRR must belong to [φ2(kmin),φ2(kmax)].

    Similar to the analysis for the scaled Robin condition, an approximate convergence factor is needed for analyzing the local maximum points. We thus introduce the following approximation of ρRR(k,p,q,v):

    ρappRR(k,p,q,v)=(aappp)2+b2app(aapp+pv)2+b2app(aappq)2+b2app(aapp+qv)2+b2app,

    where aapp and bapp are defined below (3.5).

    Theorem 4.2. (Approximation error) There exists a constant C depending on the parameters k,p, and v but independent of τ such that the following estimate holds:

    |ρRR(k,p,q,v)ρappRR(k,p,q,v)|Cτ.

    Proof. The proof is very similar to that for Theorem 3.2.

    Lemma 4.2. For small τ, the convergence factor ρRR(k,p,q,v) asymptotically achieves its maxima at kmin, kmax or ¯k=pq, i.e.,

    maxkK+ρRR(k,p,q,v)=max{ρappRR(kmin,p,q,v),ρappRR(kmax,p,q,v),ρappRR(¯k,p,q,v)}+O(τ). (4.5)

    Proof. Taking the derivative of ρappRR(k,p,q,v) in k gives

    signρappRRk=sign((kpq)g(k)), (4.6)

    where

    g(k)=2(v1)(pvq)k52+2k2v+2(q2pq(v24v+1)p2v2)k2p2q2(v1)(pvq)k+2p2q2v. (4.7)

    Investigating g(k) defined in (4.7), one finds from (4.6) that ρappRR have two interior minimum points k2>k1>0 and a unique interior maximum point ¯k=pq. Thus, ρappRR must achieve its maximum at one of the following local maximum points: kmin, kmax, and ¯k. The asymptotic result (4.5) then follows from Theorem 4.2.

    Now, we are ready to solve the min-max problem (4.2) by equi-oscillating the three maxima at kmin, kmax, and ¯k, which can be analyzed by following the idea in [32], with different details. First, we derive the equi-oscillation at the two endpoints.

    Theorem 4.3. The solution (p,q) to the min-max problem (4.2) must satisfy the equi-oscillation at kmin and kmax,

    ρRR(kmin,p,q,v)=ρRR(kmax,p,q,v). (4.8)

    Proof. It is not hard to check that for (p,q)[φ1(kmin),φ1(kmax)]×[φ2(kmin),φ2(kmax)] one has

    ρRR(kmin,p,q,v)p>0,ρRR(kmin,p,q,v)q>0,
    ρRR(kmax,p,q)p<0,ρRR(kmax,p,q)q<0.

    Thus, increasing p or q would increase ρRR(kmin,p,q,v) and decrease ρRR(kmax,p,q,v) simultaneously. While at k=¯k it holds that

    signρRR(¯k,p,q,v)p=sign(pφ1(¯k)),signρRR(¯k,p,q,v)q=sign(qφ2(¯k)),

    which is hard to be directly applied because φ1(¯k) and φ2(¯k) are different. Noting from Theorem 4.2 the fact that ρRR(¯k,p,q,v)=ρappRR(¯k,p,q,v)+O(τ), it is easy to find

    signρappRR(¯k,p,q,v)p=signρappRR(¯k,p,q,v)q=sign(2pq(pvq)2+pq(1v)).

    We thus conclude that ρRR(¯k,p,q,v) has opposite monotonicity in p and q for sufficiently small τ. Now the convergence factor ρRR can be improved until equi-oscillation is achieved as follows. Without loss of generality, we assume ρRR(kmin,p,q,v)<ρRR(kmax,p,q,v) and consider the effect of varying parameters p,q on the convergence factor ρRR. Note that the analysis starts from p<q for fixed q, and we make no assumption on ρRR(¯k,p,q,v).

    We first consider increasing p gradually. Noting that increasing p will decrease both ρRR(¯k,p,q,v) and ρRR(kmax,p,q,v), and simultaneously increase ρRR(kmin,p,q,v) (see the left panel of Figure 1), if in this process it always holds that ρRR(¯k,p,q,v)ρRR(kmax,p,q,v), then when ρRR(kmin,p,q,v)=ρRR(kmax,p,q,v) is reached, the target result is achieved, as shown in the right plot of Figure 1.

    Figure 1.  The evolution sketch of the convergence factor ρRR as increasing p for p<q.

    If the target result is not attained in the previous operation, there are two possibilities: the first is, for p<q, that ρRR(¯k,p,q,v)>ρRR(kmax,p,q,v) happens at some p before the equi-oscillation at the endpoints, and the second is that even when p increases to q (i.e., p=q), it still holds that ρRR(kmin,p,q,v)<ρRR(kmax,p,q,v) and ρRR(¯k,p,q,v)<ρRR(kmax,p,q,v). For the first possibility, we decrease q gradually such that ρRR(¯k,p,q,v)=ρRR(kmax,p,q,v) and p<q hold. Now, we return to the previous operation and continue increasing p. If ρRR(¯k,p,q,v)>ρRR(kmax,p,q,v) occurs again in the process of increasing p, we decrease q again, as mentioned above. For the second possibility, we analyze the case p=q. Noting that

    ρappRR(¯k,p,q,v)p=0,2ρappRR(¯k,p,q,v)p2=(v2+2vv+1)(v2+2v+1)3>0,

    together with Theorem 4.2, we conclude that ρRR(¯k,p,q,v) achieves its local minimum in p in an asymptotic sense. Now we continue increasing p such that p>q. In this process, increasing p will increase both ρRR(kmin,p,q,v) and ρRR(¯k,p,q,v), and simultaneously decrease ρRR(kmax,p,q,v) (see the first panel of Figure 2), if we have ρRR(kmin,p,q,v)=ρRR(kmax,p,q,v)ρRR(¯k,p,q,v) (see the third panel of Figure 2), and then the convergence factor has been improved.

    Figure 2.  The evolution sketch of the convergence factor ρRR as increasing p or q for p>q.

    If ρRR(¯k,p,q,v)=ρRR(kmax,p,q,v) is achieved for some p before the endpoints equi-oscillation, we then fix the current p and increase the other parameter q to improve the convergence factor, as sketched in the second panel of Figure 2. If it always holds that ρRR(¯k,p,q,v)ρRR(kmax,p,q,v) in this process, we arrive at our result when ρRR(kmin,p,q,v)=ρRR(kmax,p,q,v), as shown in the third plot of Figure 2.

    Otherwise if ρRR(¯k,p,q,v)>ρRR(kmax,p,q,v) in the process of increasing q for q<p, we gradually decrease p such that ρRR(¯k,p,q,v)=ρRR(kmax,p,q,v) holds. Then we still need to continue increasing q and repeat the operation mentioned above.

    We claim that in the process of continuously improving the convergence factor, ρRR(kmin,p,q,v)=ρRR(kmax,p,q,v) must be reached, since when p is approaching φ1(kmax) or q is approaching φ2(kmax), noting that both φ1(kmax) and φ2(kmax) tend to infinity for sufficiently large kmax, we must have ρRR(kmin,p,q,v)>ρRR(kmax,p,q,v). At the same time, the convergence factor ρRR is surely improved until equi-oscillation at the endpoints, and the max[kmin,kmax]ρRR is decreasing along the process.

    Theorem 4.3 shows that the solution to the min-max problem (4.2) must satisfy the endpoints equi-oscillation ρRR(kmin,p,q,v)=ρRR(kmax,p,q,v), which determines via the implicit function theorem the relation between parameters p and q, i.e., q=q(p).

    Lemma 4.3. The function q=q(p) determined by the endpoints equi-oscillation (4.8) decreases strictly in p[φ1(kmin),φ1(kmax)].

    Proof. For ease of discussion, we define F1(p,q):=ρRR(kmin,p,q,v) and F2(p,q):=ρRR(kmax,p,q,v). Then, there must exist a pair of values (p,q) such that F(p,q):=F1(p,q)F2(p,q)=0 from (4.8). Since F(p,q(p))=0, taking the derivative in p by the chain rule gives

    dF(p,q(p))dp=dF1(p,q(p))dF2(p,q(p))dp=F1F2p+F1F2qq(p)=0,

    which reveals

    q(p)=(F2pF1p)/(F1qF2q). (4.9)

    Analyzing the signs of each term gives F1p>0, F1q>0 and F2q<0, F2q<0, and we then conclude that q(p)<0 for any p[φ1(kmin),φ1(kmax)], thus q(p) is a strictly decreasing function in p.

    Lemma 4.3 shows that the two-parameter min-max problem (4.2) can be reduced to a single parameter one, which could be solved by equi-oscillation at two points, say ¯k and kmax, as shown later. In addition, we can restrict the free parameter p to a more narrow interval [q1(φ2(kmax)),φ1(kmax)], where the term q1(φ2(kmax)) is calculated from the restricted interval of q.

    Now the monotonicity of ρRR(¯k,p,q(p),v) and ρRR(kmax,p,q(p),v) in p should be investigated, which is used to solve the min-max problem

    minp[q1(φ2(kmax)),φ1(kmax)]max{ρRR(¯k,p,q(p),v),ρRR(kmax,p,q(p),v)}. (4.10)

    Lemma 4.4. At both k=kmax and k=¯k, the convergence factor ρRR(k,p,q(p),v) achieves its local maxima in p at the endpoints of [q1(φ2(kmax)),φ1(kmax)]. At k=¯k the convergence factor ρRR(k,p,q(p),v) attains its local minimum in p at φ1(¯k).

    Proof. Define F2(p,q(p)):=ρRR(kmax,p,q(p),v) as in the proof of Lemma 4.3. Analyzing the monotonicity of F2(p,q(p)) in p by directly using the derivative of F2(p,q(p)) is difficult, so we just compute the derivatives at the endpoints of [q1(φ2(kmax)),φ1(kmax)] respectively,

    dF2(q1(φ2(kmax)),φ2(kmax))dp=F2(q1(φ2(kmax)),φ2(kmax))p+F2(q1(φ2(kmax)),φ2(kmax))qq(q1(φ2(kmax)))<0,

    and

    dF2(φ1(kmax),q(φ1(kmax)))dp=F2(φ1(kmax),q(φ1(kmax)))p+F2(φ1(kmax),q(φ1(kmax)))qq(φ1(kmax))>0.

    Hence we conclude that F2(p,q(p)) is decreasing near p=q1(φ2(kmax)), while near p=φ1(kmax), F2(p,q(p)) is increasing, which shows the assertion at kmax.

    For the rest of the assertion, we now define ¯F(p,q(p)):=ρRR(¯k,p,q(p),v) and explore the monotonicity behavior of ¯F(p,q(p)) in p. Taking the derivative of ¯F(p,q(p)) in p, one finds d¯F(p,q(p))dp=0 if p=φ1(¯k) (at the same time q=q(φ1(¯k))=φ2(¯k)), which separates the problem into two cases: p>φ1(¯k) and p<φ1(¯k). When p>φ1(¯k) (i.e., q(p)<φ2(¯k)), one finds

    d¯F(p,q(p))dp=¯F(p,q(p))p+¯F(p,q(p))qq(p)>0,

    hence ¯F(p,q(p)) is strictly increasing for p(φ1(¯k),φ1(kmax)]. While for the second case p<φ1(¯k) (i.e., q(p)>φ2(¯k)), ¯F(p,q(p)) is strictly decreasing for p[q1(φ2(kmax)),φ1(¯k)) since

    d¯F(p,q(p))dp=¯F(p,q(p))p+¯F(p,q(p))qq(p)<0.

    We then arrive at our assertion at ¯k.

    Next, a simple analysis shows that, for different values of v, the best transmission parameters (p,q(p)) are located at different regions determined by p>φ1(¯k) and p<φ1(¯k), see the following lemma.

    Lemma 4.5. For v>1, the optimized parameter pair (p,q(p)) to the min-max problem (4.2) satisfies q(p)<φ1(¯k;v)<p. For v<1, the optimized parameter pair (p,q(p)) satisfies q(p)>φ1(¯k;v)>p, where we use φi(k;v) to clearly show the dependence of φi on v.

    Proof. By the definitions (4.3) and (4.4) of φ1 and φ2, it is easy to conclude that φ1(k;v)=φ2(k;1v). Thus, when v>1 it holds that φ2(¯k;v)<φ2(¯k;1v)=φ1(¯k;v), since φ2(k;v) is decreasing in v:

    φ2(k;v)v=a2+b22a(a2(v+1)+b2(v1)(a2+b2)(a2(v+1)2+b2(v1)2)1)<0.

    In addition, we recall the expression of the convergence factor

    ρRR(k,p,q(p),v)=(ap)2+b2(a+pv)2+b2(aq(p))2+b2(a+q(p)v)2+b2

    and noting the assumption v>1, we find that by swapping the parameters p and q(p), we get a larger convergence factor

    ρRR(k,p,q(p),v)<ρRR(k,q(p),p,v).

    Thus, for v>1 the optimizer can only be obtained for p>q(p). To compare p and φ1(¯k;v), noting the proved fact that increasing p would increase ρRR(¯k,p,q(p),v) for p>q(p), and the montonicity of ρRR(¯k,p,q(p),v) in p is determined by the sign of pφ1(¯k;v), we infer from the proof of the case p>q in Theorem 4.3 that p>φ1(¯k;v). This proves the first assertion and q(p)<φ2(¯k;v)<φ1(¯k;v)<p holds. The second assertion for v<1 can be obtained in a similar fashion.

    Remark 4.1. Lemma 4.5 shows, for the min-max problem (4.10), that we are urged to consider only the case v>1, and the other case v<1 reduces to the former one by swapping the value of p and q(p) and replacing v by 1v. For v>1, the min-max problem (4.10) is equivalently reduced to

    minφ1(¯k)<pφ1(kmax)max{¯F(p,q(p)),F2(p,q(p))}. (4.11)

    Theorem 4.4. (Solve by equi-oscillation) For v>1 and kmax large enough, the parameter pair (p,q) with p>q determined by equi-oscillating the convergence factor ρRR at ¯k and kmax, i.e.,

    ρRR(¯k,p,q(p),v)=ρRR(kmax,p,q(p),v), (4.12)

    is the unique solution to the min-max problem (4.11).

    Proof. Indicated by Lemma 4.5, we only need to consider the case p>φ1(¯k;v)>q(p). To solve the min-max problem (4.11), we compare the values of ¯F(p,q(p)) and F2(p,q(p)) (the convergence factor evaluated at ¯k and kmax) at the endpoints of the searching interval [φ1(¯k),φ1(kmax)] of p. Before proceeding with this, we need to expand φ1(kmax) for sufficiently large kmax,

    φ1(kmax)=5v28v+54v(v1)τ+2(v1)τ32vk2max+O(1k2max).

    Based on this result, we find the values at the right boundary φ1(kmax) of p indicated by Lemma 4.4:

    ¯F(φ1(kmax),q(φ1(kmax)))=q(φ1(kmax))2τ+1q(φ1(kmax))2τ+v2+O(1k2max),F2(φ1(kmax),q(φ1(kmax)))=1v+O(1k2max),

    which show that F2(φ1(kmax),q(φ1(kmax)))<¯F(φ1(kmax),q(φ1(kmax))) for kmax large enough. In addition, it is easy to check that F2(φ1(¯k),φ2(¯k))>¯F(φ1(¯k),φ2(¯k)). Hence, the solution (p,q) to the min-max problem (4.11) is determined by the minimum of the intersections between ¯F(p,q(p)) and F2(p,q(p)). Now, to end the proof, we only need to show that the intersection of the above two functions is unique. To do this, we define L(p):=¯F(p,q(p))F2(p,q(p)) and recall q(p)<0 by (4.9). Direct calculations show that ¯FpF2p>0 and ¯FqF2q<0 for any p>φ1(¯k). These findings indicate that

    dLdp=¯FpF2p+(¯FqF2q)q(p)>0.

    We then conclude that L(p) is strictly increasing in p for p>φ1(¯k), starting from a negative value and strictly increasing to a positive value. The uniqueness of the interaction between ¯F(p,q(p)) and F2(p,q(p)) has arrived and ¯F(p,q(p))=F2(p,q(p)) defines the unique solution to the min-max problem (4.11).

    Theorem 4.5. For v>1, the min-max problem (4.1) is solved in an asymptotic sense of kmax large by

    p=2(v1)2vk12max+234a12min(v2+1)v(v1)k14max,q=2aminvv1254v((a2min+b2min)(v2+1)+2v(a2minb2min))amin(v1)3k14max, (4.13)

    and the corresponding convergence factor ρRR satisfies the following estimate:

    maxkminkkmaxρ(k,p,q,v)=1v254(v+1)aminv(v1)k14max+O(k12max). (4.14)

    Proof. Based on Theorems 4.1, 4.3, and 4.4, Lemma 4.3, and Remark 4.1, we only need to show that the optimized parameter pair (p,q) solves the equi-oscillation problem (4.12). For v>1, guided by numerical investigations, we make the ansatz p=c1k14max+c2k12max and q=c3c4k14max. Using Lemma 4.2 one then finds ¯k=pq=c2c3k12max+(c1c3c2c4)k14maxc1c4. Applying these ansatz in ρRR(kmin,p,q,v), ρRR(kmax,p,q,v), and ρRR(¯k,p,q,v) and expanding for large kmax gives

    ρRR(kmin,p,q,v)=(c3amin)2+b2min(c3+aminv)2+b2minv2(a3minv+c3a2min(v1)(b2minvc23)amin+c3b2min(v1))(v+1)c4(c3amin)2+b2min((aminv+c3)2+b2minv2)32k14max+O(k12max), (4.15)
    ρRR(kmax,p,q,v)=c222c2+1c22v2+2c2v+1c1(v+1)(2c2v2c22c22v+2)2c222c2+1(c22v2+2c2v+1)32k14max+O(k12max) (4.16)

    and

    ρRR(¯k,p,q,v)=1v2c3(v+1)c2v2k14max+O(k12max). (4.17)

    Now, equalizing (4.15)–(4.17), that is, matching the first two terms, noting the conditions v>1 and c1,c2,c3,c4>0, we arrive at our result (4.13), and (4.14) follows.

    Remark 4.2. The convergence factor estimate (4.14) shows that the OSWR algorithm with the Robin-Robin transmission condition achieves mesh-independent convergence if the mesh size is not too large. This is a property that cannot be reached by the scaled Robin condition, since the requirement of p=q in the scaled Robin condition deviates significantly from the optimum (4.13). Additionally, similar to Remark 3.1, one can conclude that a large heterogeneity contrast v leads to faster convergence of the OSWR algorithm. Besides, the Robin-Robin transmission also applies to homogeneous media β1=β2, however, loses the mesh-independent convergence. In fact, the results in Theorem 4.5 cannot be directly applied since p and q are singular at v=1. Using a similar but much simpler analysis, the optimized parameters for v=1 are obtained as

    p=258a14mink38max,q=218a34mink18max

    and the corresponding convergence factor estimate is

    maxkminkkmaxρRR(k,p,q,v)=1298a14mink18max+O(k14max).

    Remark 4.3. Our analysis of the Robin-Robin condition simplifies to the standard Fourier heat conduction in composites when τ=0, as shown in Remark 3.2 for the scaled Robin condition. Additionally, the introduction of the phase-lag τ also slows down the convergence of the OSWR algorithm, as indicated by the convergence factor estimate (4.14). However, in contrast to the significant slowdown in convergence caused by a relatively large phase-lag when using the scaled Robin condition, the phase-lag only has a very small effect on the convergence of the OSWR algorithm when the Robin-Robin condition is applied. This is due to the mesh-independent convergence property of the Robin-Robin condition, which relies solely on the heterogeneity contrast between different media.

    In this section, to illustrate our theoretical results, we perform several numerical experiments for our model problem (1.2) on the space-time domain Ω×(0,T]=[1,1]×(0,1], with nonhomogeneous Dirichlet conditions at x=1 and x=1. The domain Ω is decomposed into two nonoverlapping subdomains Ω=Ω1Ω2 with Ω1=[1,0) and Ω2=(0,1] according to heterogeneous media. If not specified, we consider only the case where the heterogeneous parameter v>1(β2>β1). To discretize the model problem in each subdomain, we use the Crank-Nicolson difference scheme for temporal discretization, and the centered finite difference method in space. The model problem is solved using the OSWR algorithm, starting from a random initial guess on the interface {0}×(0,T] (such that all the possible frequency components are present) and stops if

    max{gn+11gn1,gn+12gn2}<106.

    Using the numerical experiments, we would like to justify the following concerns.

    1) Does the asymptotic formula (3.12) ((4.13), resp.) well predict the optimal parameters to be applied in the scaled Robin (Robin-Robin, resp.) condition?

    2) The Robin-Robin condition leads the OSWR algorithm to a mesh-independent convergence behavior (see Theorem 4.5), while the scaled Robin can only lead to mesh-dependent convergence (Theorem 3.4).

    3) Larger heterogeneity contrast leads to faster convergence, for both transmission conditions, as indicated by Theorems 3.4 and 4.5.

    4) The phase-lag time τ deteriorates the algorithm's performance: the larger the phase-lag τ is, the slower the OSWR algorithm converges, as indicated by Remarks 3.2 and 4.3. However, the Robin-Robin condition is more robust to the phase-lag time than the scaled Robin condition.

    To examine the first concern, let τ1=τ2=0.001, β1=0.0001, and β2=0.5. We consider our model problem (1.2) with the exact solution [23]

    T(x,t)={154x2(x21)t2etifx0,4x2sin(3t4)ifx>0, (5.1)

    which is also used to determine the initial functions T0,˜T0 and the source term f. To check the usefulness of our prediction, for the scaled Robin transmission condition, we sample the parameter p near the optimized parameter p in (3.12) and count the number of iterations required by the OSWR algorithm, with a similar process for the Robin-Robin transmission condition. We then plot the required number of iterations in Figure 3. The results show that for both transmission conditions, our asymptotic predictions are very close to the numerical optimal. The difference may be attributed to the application of Fourier analysis in time, see [46] for more details.

    Figure 3.  Number of iterations corresponds to the optimized parameters (red ) compared with other parameter values. The scaled Robin condition is on the left, and the Robin-Robin condition is on the right. A uniform mesh of size Δx=Δt=0.002 is applied.

    To justify the second concern, we still consider the model problem (1.2) with the exact solution (5.1). The phase-lag time is again set as τ1=τ2=0.001, with the parameters β1=0.0001, β2=0.04. For different time step sizes Δt, we record the number of iterations required by the OSWR algorithms (2.1) and (2.2) and present them in Table 1. From these results, one easily finds that the Robin-Robin condition converges much faster than the scaled Robin condition. More importantly, the number of iterations required by the Robin-Robin condition does not vary as the mesh size decreases, implying the mesh-independent convergence behavior as predicted by Theorem 4.5.

    Table 1.  Number of iterations required by the OSWR algorithms (2.1) and (2.2) to reach a tolerance of 106 for different mesh size Δt.
    Δt 132 164 1128 1256 1512
    Scaled Robin 11 13 13 15 17
    Robin-Robin 9 9 9 11 11

     | Show Table
    DownLoad: CSV

    We now look at how the heterogeneity contrast ν influences the convergence behavior. To this end, we still consider the exact solution mentioned above and set τ1=τ2=0.001. Now, we fix β1=0.0001 and vary β2 such that the heterogeneity contrast ν varies. In Figure 4 we plot the relative error reduction as a function of the number of iterations, where we see that for both the scaled Robin and Robin-Robin transmission conditions, larger heterogeneity contrast requires a smaller number of iterations, as shown in Theorems 3.4 and 4.5.

    Figure 4.  For different heterogeneity contrast ν, the relative error reduction is shown as a function of the number of iterations for the OSWR algorithms (2.1) and (2.2) when using a uniform mesh of size Δx=Δt=0.002. Left: the scaled Robin condition; Right: the Robin-Robin condition.

    To illustrate the last concern, we still consider the exact solution (5.1) with β1=0.002 and β2=0.5. Table 2 and Figure 5 show the number of iterations required by the OSWR algorithms (2.1) and (2.2) for different phase-lag time τ, which indicates that for both transmission conditions a larger phase-lag time would lead to a deterioration of the convergence: the larger the phase-lag time, the slower the algorithm converges, as indicated by Remarks 3.2 and 4.3. In addition, the Robin-Robin condition outperforms the scaled Robin for all phase-lag τ, especially when the phase-lag is relatively large, showing its robustness to the phase-lag over the scaled Robin condition. However, we do not observe a phase-lag independent performance for the Robin-Robin condition. This is because, in our opinion, the mesh size used is not small enough such that the asymptotic regime in (4.14) is reached.

    Table 2.  Number of iterations required by the OSWR algorithms (2.1) and (2.2) is shown using scaled Robin and Robin-Robinconditions with a uniform mesh of size Δx=Δt=0.002, for different phase-lag times τ.
    τ Scaled Robin Robin-Robin
    0 17 11
    0.0002 19 13
    0.0004 23 15
    0.0008 31 19
    0.0016 50 24
    0.0032 130 35

     | Show Table
    DownLoad: CSV
    Figure 5.  The plot of the number of iterations shown in the Table 2 as a function of the phase-lag τ.

    In addition, though our analysis is based on the assumption τ1=τ2, the theoretical results still hold for heterogeneous τ. To illustrate this result, we still consider the exact solution (5.1) with β1=0.0001 and β2=0.04. Setting S1=β2p(τ1)I and S2=β1p(τ2)I in the scaled Robin and S1=β2p(τ1)I and S2=β1q(τ2)I in the Robin-Robin conditions, we collect the number of iterations required by the OSWR algorithms (2.1) and (2.2) in Table 3, which shows the effectiveness of our algorithm in solving problems with heterogeneous phase-lag τ.

    Table 3.  Number of iterations required by the OSWR algorithms (2.1) and (2.2) to reach a tolerance of 106 for different mesh size Δt when τ1=0.001 and τ2=0.0005.
    Δt 132 164 1128 1256 1512
    Scaled Robin 11 13 13 15 16
    Robin-Robin 9 9 9 11 11

     | Show Table
    DownLoad: CSV

    We propose to solve the heterogeneous non-Fourier heat transfer problem using the OSWR method, where two types of Robin transmission conditions are proposed because the phase-lag time τ is generally small. We have rigorously optimized the convergence factors obtained using Fourier analysis for both transmission conditions through asymptotic analysis. The results reveal several interesting new findings: for both conditions, the large heterogeneity contrast benefits the convergence; notably, the OSWR with the Robin-Robin condition is asymptotically mesh-independent, while the OSWR with the scaled Robin condition is not; more importantly, compared to standard Fourier heat transfer, the phase-lag time slows down the convergence speed, especially when the phase-lag time is relatively large, which poses a challenging topic for further investigation, particularly when the mesh size is relatively large.

    Feng Hu: Conceptualization, Formal analysis, Investigation, Methodology, Software, Visualization, Validation, Writing-original draft and editing; Yingxiang Xu: Conceptualization, Formal analysis, Investigation, Methodology, Supervision, Project administration, Validation, and Writing-review and editing.

    The authors declare they have used Artificial Intelligence (AI) tools in the creation of this article.

    This work is supported in part by the Science and Technology Development Planning of Jilin Province under grant YDZJ202201ZYTS573, the National Key R & D Program of China under grant 2021YFA1003400, the National Natural Science Foundation of China under grant 12071069, and the Fundamental Research Funds for the Central Universities under grant 2412022ZD032.

    There is no known conflict of interest associated with this publication.



    [1] A. Narasimhan, S. Sadasivam, Non-Fourier bio heat transfer modelling of thermal damage during retinal laser irradiation, Int. J. Heat Mass Tran., 60 (2013), 591–597. https://doi.org/10.1016/j.ijheatmasstransfer.2013.01.010 doi: 10.1016/j.ijheatmasstransfer.2013.01.010
    [2] J. Zhou, Y. Zhang, J. Chen, Non-Fourier heat conduction effect on laser-induced thermal damage in biological tissues, Numer. Heat Tr. A-Appl., 54 (2008), 1–19. https://doi.org/10.1080/10407780802025911 doi: 10.1080/10407780802025911
    [3] T. Xue, X. Zhang, K. K. Tamma, Investigation of thermal interfacial problems involving non-locality in space and time, Int. Commun. Heat Mass, 99 (2018), 37–42. https://doi.org/10.1016/j.icheatmasstransfer.2018.10.008 doi: 10.1016/j.icheatmasstransfer.2018.10.008
    [4] Z. Y. Guo, Y. S. Xu, Non-Fourier heat conduction in IC chip, J. Electron. Packag. Sep., 117 (1995), 174–177. https://doi.org/10.1115/1.2792088 doi: 10.1115/1.2792088
    [5] H. D. Wang, B. Y. Cao, Z. Y. Guo, Non-Fourier heat conduction in carbon nanotubes, J. Heat Transfer, 134 (2012), 051004. https://doi.org/10.1115/1.4005634 doi: 10.1115/1.4005634
    [6] C. Cattaneo, A form of heat-conduction equations which eliminates the paradox of instantaneous propagation, Comptes Rendus, 247 (1958), 431.
    [7] P. Vernotte, Some possible complications in the phenomena of thermal conduction, Compte Rendus, 252 (1961), 2190–2191.
    [8] M. Ozisik, D. Y. Tzou, On the wave theory in heat conduction, J. Heat Transfer, 116 (1994), 526–535. https://doi.org/10.1115/1.2910903 doi: 10.1115/1.2910903
    [9] J. Ordonez-Miranda, J. Alvarado-Gil, Thermal wave oscillations and thermal relaxation time determination in a hyperbolic heat transport model, Int. J. Therm. Sci., 48 (2009), 2053–2062. https://doi.org/10.1016/j.ijthermalsci.2009.03.008 doi: 10.1016/j.ijthermalsci.2009.03.008
    [10] A. Vedavarz, S. Kumar, M. K. Moallemi, Significance of non-Fourier heat waves in conduction, J. Heat Transfer, 116 (1994), 221–224. https://doi.org/10.1115/1.2910859 doi: 10.1115/1.2910859
    [11] D. Y. Tzou, J. C. Dowell, Computational techniques for microscale heat transfer, In: Handbook of Numerical Heat Transfer, 2000. https://doi.org/10.1002/9780470172599.ch20
    [12] K. C. Liu, P. J. Cheng, Finite propagation of heat transfer in a multilayer tissue, J. Thermophys. Heat Tr., 22 (2008), 775–782. https://doi.org/10.2514/1.37267 doi: 10.2514/1.37267
    [13] C. Barman, P. Rath, A. Bhattacharya, A non-Fourier bioheat transfer model for cryosurgery of tumor tissue with minimum collateral damage, Comput. Meth. Prog. Bio., 200 (2021), 105857. https://doi.org/10.1016/j.cmpb.2020.105857 doi: 10.1016/j.cmpb.2020.105857
    [14] P. Namakshenas, A. Mojra, Microstructure-based non-Fourier heat transfer modeling of HIFU treatment for thyroid cancer, Comput. Meth. Prog. Bio., 197 (2020), 105698. https://doi.org/10.1016/j.cmpb.2020.105698 doi: 10.1016/j.cmpb.2020.105698
    [15] Z. S. Deng, J. Liu, Non-Fourier heat conduction effect on prediction of temperature transients and thermal stress in skin cryopreservation, J. Therm. Stresses, 26 (2003), 779–798. https://doi.org/10.1080/01495730390219377 doi: 10.1080/01495730390219377
    [16] A. K. Kheibari, M. Jafari, M. B. Nazari, Propagation of heat wave in composite cylinder using Cattaneo-Vernotte theory, Int. J. Heat Mass Transfer, 160 (2020), 120208. https://doi.org/10.1016/j.ijheatmasstransfer.2020.120208 doi: 10.1016/j.ijheatmasstransfer.2020.120208
    [17] A. Fehér, R. Kovacs, Novel evaluation method for non-Fourier effects in heat pulse experiments, 2021. https://doi.org/10.48550/arXiv.2101.01123
    [18] A. Fehér, D. Markovics, T. Fodor, R. Kovacs, Size effects and non-Fourier thermal behaviour in rocks, ISRM EUROCK, 2020.
    [19] P. Duhamel, A new finite integral transform pair for hyperbolic conduction problems in heterogeneous media, Int. J. Heat Mass Transfer, 44 (2001), 3307-3320. https://doi.org/10.1016/S0017-9310(00)00360-4 doi: 10.1016/S0017-9310(00)00360-4
    [20] S. Singh, Z. Li, A high order compact scheme for a thermal wave model of bio-heat transfer with an interface, Numer. Math. Theory Me., 11 (2018), 321–337. https://doi.org/10.4208/nmtma.OA-2017-0048 doi: 10.4208/nmtma.OA-2017-0048
    [21] M. Asif, F. Bilal, R. Bilal, N. Haider, S. A. M. Abdelmohsenc, S. M. Eldind, An efficient algorithm for the numerical solution of telegraph interface model with discontinuous coefficients via Haar wavelets, Alex. Eng. J., 72 (2023), 275–285. https://doi.org/10.1016/j.aej.2023.03.074 doi: 10.1016/j.aej.2023.03.074
    [22] T. Ahmod, J. Dutta, Finite element method for hyperbolic heat conduction model with discontinuous coefficients in one dimension, Proc. Math. Sci., 132 (2022), 6. https://doi.org/10.1007/s12044-021-00646-3 doi: 10.1007/s12044-021-00646-3
    [23] B. Deka, J. Dutta, Finite element methods for non-Fourier thermal wave model of bio heat transfer with an interface, J. Appl. Math. Comput., 62 (2020), 701–724. https://doi.org/10.1007/s12190-019-01304-8 doi: 10.1007/s12190-019-01304-8
    [24] B. Deka, J. Dutta, Convergence of finite element methods for hyperbolic heat conduction model with an interface, Comput. Math. Appl., 79 (2020), 3139–3159. https://doi.org/10.1016/j.camwa.2020.01.013 doi: 10.1016/j.camwa.2020.01.013
    [25] J. Dutta, B. Deka, Optimal a priori error estimates for the finite element approximation of dual-phase-lag bio heat model in heterogeneous medium, J. Sci. Comput., 87 (2021), 58. https://doi.org/10.1007/s10915-021-01460-9 doi: 10.1007/s10915-021-01460-9
    [26] S. Han, Finite volume solution of 2-D hyperbolic conduction in a heterogeneous medium, Numer. Heat Tr. A-Appl., 70 (2016), 723–737. https://doi.org/10.1080/10407782.2016.1193347 doi: 10.1080/10407782.2016.1193347
    [27] H. Sauerland, T. P. Fries, The stable XFEM for two-phase flows, Comput. Fluids, 87 (2013), 41–49. https://doi.org/10.1016/j.compfluid.2012.10.017 doi: 10.1016/j.compfluid.2012.10.017
    [28] B. Ayuso de Dios, M. Holst, Y. Zhu, L. Zikatanov, Multilevel preconditioners for discontinuous Galerkin approximations of elliptic problems with jump coefficients, Math. Comp., 83 (2014), 1083–1120. https://doi.org/10.1090/S0025-5718-2013-02760-3 doi: 10.1090/S0025-5718-2013-02760-3
    [29] M. J. Gander, Optimized Schwarz methods, SIAM J. Numer. Anal., 44 (2006), 699–731. https://doi.org/10.1137/S0036142903425409 doi: 10.1137/S0036142903425409
    [30] Y. Maday, F. Magoules, Optimized Schwarz methods without overlap for highly heterogeneous media, Comput. Method. Appl. M., 196 (2007), 1541–1553. https://doi.org/10.1016/j.cma.2005.05.059 doi: 10.1016/j.cma.2005.05.059
    [31] M. J. Gander, O. Dnbois, Optimized Schwarz methods for a diffusion problem with discontinuous coefficient, Numer. Algor., 69 (2015), 109–144. https://doi.org/10.1007/s11075-014-9884-2 doi: 10.1007/s11075-014-9884-2
    [32] M. J. Gander, T. Vanzan, Heterogeneous optimized Schwarz methods for second order elliptic PDEs, SIAM J. Sci. Comput., 41 (2019), A2329–A2354. https://doi.org/10.1137/18M122114X doi: 10.1137/18M122114X
    [33] M. J. Gander, S. B. Lunowa, Ch. Rohde, Non-overlapping Schwarz Waveform-Relaxation for Nonlinear Advection-Diffusion Equations, SIAM J. Sci. Comput., 45 (2023), A49–A73. https://doi.org/10.1137/21M1415005 doi: 10.1137/21M1415005
    [34] M. J. Gander, L. Halpern, M. Kern, A Schwarz waveform relaxation method for advection-diffusion-reaction problems with discontinuous coefficients and non-matching grids, In: Domain Decomposition Methods in Science and Engineering, Berlin: Springer, 2007,283–290. https://doi.org/10.1007/978-3-540-34469-8_33
    [35] L. Halpern, C. Japhet, Discontinuous Galerkin and nonconforming in time optimized Schwarz waveform relaxation for heterogeneous problems, In: Domain Decomposition Methods in Science and Engineering XVII, Berlin: Springer, 2008,211–219. https://doi.org/10.1007/978-3-540-75199-1_23
    [36] L. Halpern, C. Japhet, J. Szeftel, Optimized Schwarz waveform relaxation and discontinuous Galerkin time stepping for heterogeneous problems, SIAM J. Numer. Anal., 50 (2012), 2588–2611. https://doi.org/10.1137/120865033 doi: 10.1137/120865033
    [37] M. D. Al-Khaleel, M. J. Gander, and P. M. Kumbhar, Optimized Schwarz waveform relaxation methods for the telegrapher equation, SIAM J. Sci. Comput., 46 (2024), A3528–A3551. https://doi.org/10.1137/24M1642962 doi: 10.1137/24M1642962
    [38] M. J. Gander, J. Lin, S. L. Wu, X. Yue, T. Zhou, Parareal: parallel-in-time algorithms based on the diagonalization technique, 2020. https://doi.org/10.48550/arXiv.2005.09158
    [39] M. J. Gander, S. Vandewalle, Analysis of the Parareal time-parallel time integration method, SIAM J. Sci. Comput., 29 (2007), 556–578. https://doi.org/10.1137/05064607X doi: 10.1137/05064607X
    [40] G. Ciaramella, M. J. Gander, I. Mazzieri, Unmapped tent pitching schemes by waveform relaxation, In: Decomposition Methods in Science and Engineering XXVII, Cham: Springer, 2024,455–462. https://doi.org/10.1007/978-3-031-50769-4_54
    [41] M. J. Gander, L. Halpern, Optimized Schwarz waveform relaxation methods for advection-reaction-diffusion problems, SIAM J. Numer. Anal., 45 (2007), 666–697. https://doi.org/10.1137/050642137 doi: 10.1137/050642137
    [42] D. Bennequin, M. J. Gander, L. Gouarin, L. Halpern, Optimized Schwarz waveform relaxation for advection-reaction-diffusion equations in two dimensions, Numer. Math., 134 (2016), 513–567. https://doi.org/10.1007/s00211-015-0784-8 doi: 10.1007/s00211-015-0784-8
    [43] V. Martin, An optimized Schwarz waveform relaxation method for the unsteady convection-diffusion equation in two dimensions, Appl. Numer. Math., 52 (2005), 401–428. https://doi.org/10.1016/j.apnum.2004.08.022 doi: 10.1016/j.apnum.2004.08.022
    [44] Y. Xu, The influence of domain truncation on the performance of optimized Schwarz methods, Electron. T. Numer. Ana., 49 (2018), 182–209. https://doi.org/10.1553/etna_vol49s182 doi: 10.1553/etna_vol49s182
    [45] W. B. Dong, H. S. Tang, Convergence analysis of Schwarz waveform relaxation method to compute coupled advection-diffusion-reaction equations, Math. Comput. Simul., 218 (2024), 462–481. https://doi.org/10.1016/j.matcom.2023.11.026 doi: 10.1016/j.matcom.2023.11.026
    [46] M. J. Gander, V. Martin, Why Fourier mode analysis in time is different when studying Schwarz waveform relaxation, J. Comput. Phys., 491 (2023), 112316. https://doi.org/10.1016/j.jcp.2023.112316 doi: 10.1016/j.jcp.2023.112316
  • Reader Comments
  • © 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(362) PDF downloads(23) Cited by(0)

Figures and Tables

Figures(5)  /  Tables(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog