An iterative procedure to solve a coupled two-fluids turbulence model

ESAIM: M2AN 44 (2010) 693–713 ESAIM: Mathematical Modelling and Numerical Analysis DOI: 10.1051/m2an/2010015 www.esaim-m2an.org AN ITERATIVE PROCEDURE TO SOLVE A COUPLED TWO-FLUIDS TURBULENCE MODEL ∗ Tomas Chac´ on Rebollo1,St ´ ephane Del Pino2and Driss Yakoubi3 Abstract. This paper introduces a scheme for the numerical approximation of a model for two turbulent flows with coupling at an interface. We consider the variational formulation of the coupled model, where the turbulent kinetic energy equation is formulated by transposition. We prove the convergence of the approximation to this formulation for 3D flows for large turbulent viscosities and smooth enough flows, whenever bounded in W1,p Sobolev norms for plarge enough. Under the same assumptions, we show that the limit is a solution of the initial problem. Finally, we give some numerical experiments to enlighten the theoretical work. Mathematics Subject Classification. 63N30, 76M10. Received March 29, 2008. Revised August 1st, 2009. Published online February 23, 2010. 1. Introduction In this contribution we focus our attention on the modelling of the surface layer between the atmosphere and the ocean. We are interested in designing effective procedures to solve the following coupled model: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ −∇·(αi(ki)∇ui)+grad pi=fiin Ωi, ∇·ui=0inΩ i, −∇·(γi(ki)∇ki)=αi(ki)|∇ui|2in Ωi, ui=0on Γi, ki=0onΓ i, αi(ki)∂niui−pini+κ(ui−uj)|ui−uj|=0on Γ,1≤i=j≤2, ki=λ|u1−u2|2on Γ. (1.1) Whereeachtriple(ui,p i,k i) is defined in the domain Ωi,1≥i≥2. The generic point in R2,resp. inR3,is denoted by x=(x, z), resp. x=(x, y, z). Keywords and phrases. Ocean-atmosphere coupling, turbulent flows, convergence analysis, iterative method, spectral method. ∗T. Chac´on Rebollo was partially funded by EU Marie-Curie Fellowship Programme, by Spanish Government Grant MTM2006- 01275 and by Junta de Andaluc´ıa Grant P07-FQM-02538. 1Departamento de Ecuaciones Diferenciales y Analisis Numerico, Universidad de Sevilla, Spain. 2CEA, DAM, DIF, 91297 Arpajon, France. 3Laboratoire Jacques-Louis Lions, Universit´e Pierre et Marie Curie, 4 place Jussieu, 75005 Paris Cedex, France. Article published by EDP Sciences c EDP Sciences, SMAI 2010 694 T. CHAC ´ ON REBOLLO ET AL. System (1.1) is a simplified model for two stationary turbulent flows in adjacent domains, coupled by boundary conditions on the interface, such as the system atmosphere-ocean. Indeed, it is a simplified mathematical formulation of the RANS (Reynolds Averaged Navier-Stokes) model of order 1 used to simulate a stationary mean flow when convection is neglected. This kind of modelling is often used in engineering or geophysics, see for instance Bernardi et al. [6], Launder and Spalding [16], Mohammadi and Pironneau [21], Piquet [22], Wilcox [23]. In what follows, Ωi(i=1,2) are bounded domains of Rd,d=2,3, which are either convex or of class C1,1, with boundaries ∂Ωi=Γ i∪Γ,Γ=Ω1∩Ω1being the interface between the two fluids. Γ is assumed to be flat. Indeed, we assume that the so-called “rigid lid hypothesis” (introduced by Bryan in [9]) holds, an hypothesis which is standard in geophysics and oceanography. Each of the two turbulent fluids is modeled by a simplified one-equation turbulence model whose unknowns are the velocity uiand the turbulent kinetic energy (TKE) ki. In the first equation we model the generation of eddy viscosity in flow iby the term −αi(ki)∇ui. The (positive) quantity αi(ki) is the eddy viscosity. This is a simplification of the usual modelling of Reynolds Stress Tensor by Ri≃−αi(ki)∇ui+∇tui. We prefer the first expression for simplicity of mathematical analysis, although our analysis still holds for the second one. We also neglect transport effects, we intend to analyze them in a forthcoming paper. The fluids are assumed to be incompressible (second equation). In the third equation we model the generation of TKE by means of a production source term αi(ki)|∇ui|2, although the physical one should be αi(ki)|∇ui+∇tui|2. Again, we prefer the first expression for simplicity of mathematical treatment. Also the turbulent diffusion of TKE is the function γi(ki). We neglect the viscous dissipation effects, to avoid to manage an additional statistic of the turbulence (a mixing length or the turbulent dissipation ε, for instance). We assume nonslipping boundary conditions in the boundary parts Γifor simplicity (fourth and fifth equations). These in practice are replaced by wall-laws to simulate the generation of turbulence on solid boundaries. The sixth equation globally models the interaction of the two boundary layers on one and another side of the interface Γ as friction effects, by means of a set of boundary conditions similar to Manning’s law. Finally, the last equation models the production of TKE in the interface. The coefficients κiand λare positive. We assume that the turbulent diffusions αiand γibelong to W1,∞(R)andverifyαi≥ν, γi≥ν,forsome ν>0. The eddy diffusions usually are unbounded functions of the TKE of the form a+b√k,asweusein the numerical simulations reported in Section 5(see for instance [8,17,21]). But this renders the analysis much more complex even for a one-fluid turbulence model (see [18]). So we consider a simplified model, that still includes several realistic non-linear interactions. System (1.1) was studied in [3] where existence and uniqueness of small smooth solutions were proved. Spectral and Finite Element discretizations were studied in subsequent papers by the same authors and coworkers (see [4,5]). In these papers, the ability of these discretization techniques to approach the solution of model (1.1) was proved. However, in both cases the discretizations achieved consisted in fully non-linear sets of algebraic equations. Our purpose here is to derive iterative procedures to solve system (1.1) that decouples the interaction of the problem, leading to mildly non-linear problems. Let us introduce the function spaces Xi=vi∈H1(Ωi); vi=0on Γi, L2 0(Ωi)=qi∈L2(Ωi); Ωi qi=0 .(1.2) AN ITERATIVE PROCEDURE TO SOLVE A COUPLED TWO-FLUIDS TURBULENCE MODEL 695 Consider also two conjugate positive real numbers rand ri.e. 1 r+1 r=1,such that r>d. We introduce the following iterative procedure to solve (1.1): once known un i∈Xi,pn i∈L2 0(Ωi), ki∈W1,r(Ωi), i=1,2, solve: Problem 1. Obtain un+1 i∈Xi,i=1,2, such that Ωi αi(kn i)∇un+1 i:∇vidx−Ωi (∇·vi)pn+1 idx+κiΓ|un+1 i−un+1 j|(un+1 i−un+1 j)·vidτ=Ωi fi·vidτ, ∀vi∈Xi,and Ωi (∇·un+1 i)qidx=0,∀qi∈L2 0(Ωi), and Problem 2. Obtain kn+1 i∈W1,r(Ωi),i=1,2, such that kn+1 i=0 onΓ i,k n+1 i=λ|un+1 1−un+1 2|2on Γ, Ωi γi(kn i)kn+1 iϕidx=Ωi αi(kn i)|∇un+1 i|2ϕidx,∀ϕi∈W1,r 0(Ωi). Remark 1.1. We take r>dto give a sense to the equation for the kiin model (1.1). Indeed, the term αi(kn i)|∇un+1 i|2belongs to L1(Ωi), and it follows from the Sobolev Imbedding Theorem that the test function ϕi in Problem 2 belongs to L∞(Ωi), so that the right-hand member is well defined. Observe that Problem 1 is in reality non-linear due to the presence of the Manning-like source term. This is a mild non-linearity due to the monotonic nature of this term, that may be made explicit in practice if mass-lumping techniques are used. Our main result states that if the sequences (un i)nand (kn i)nare respectively bounded in W1,3+ε(Ωi)dand W1,3(Ωi), then, for small enough data (in a convenient sense), the iterative scheme is contracting. This regularity is realistic, as it is not far from the W1,2regularity that has been proved for problem (1.1) for general data. The main ingredients to show the convergence of our scheme are the convenient choices of test functions, and the use of the harmonic liftings Riof Dirichlet boundary conditions on Γ on the Ωi(see the proof of Thm. 3.4). Our paper is organised as follows. In Section 2we introduce a weak formulation of the above iterative procedure. Section 3is devoted to prove the contractiveness of the TKE sequence. Due to the production term of the TKE on interface Γ: ki=λ|u1−u2|2, it is necessary to estimate the expression |un+1 1−un+1 2|2−|un 1−un 2|2H 1 2 00(Γ) ,(1.3) where the special space H 1 2 00(Γ) is the subspace of H1 2(Γ) whose extension by zero to ∂Ω1(for instance, it could be also to ∂Ω2) belongs to H1 2(∂Ω1). An intrinsic scalar product on H 1 2 00(Γ) is defined as ((u, v))H 1 2 00 (Γ) =Γ u(x)v(x)dx+ΓΓ (u(x)−u(y)) (v(x)−v(y)) |x−y|ddxdy+Γ u(x)v(x) d(x, ∂Γ) ,(1.4) where the first two summands define the H1 2(Γ) scalar product (see Adams and Fournier [1], Thm. 7.48). Its expression involves the distance d(x, ∂Γ) to the boundary of ∂Γ. It comes from the restriction to H 1 2 00(Γ) of the scalar product in H1 2(∂Ω1) for instance. It is given by Lions and Magenes in [20], Chapter 1, Theorem 11.7. 696 T. CHAC ´ ON REBOLLO ET AL. The estimation of (1.3) is done using a Grisvard’s result (Lem. 3.5 in the paper at hand, see [14]forthe original reference) on estimates in Ws,p of products of functions of Wsj,pj, for some real numbers s, sjand some non-negative integers p, pj,j=1,2. Concerning the sequence of pressure iterates (pn i)n, we use a specific inf-sup condition (see Cor. 3.7) to show that it is a Cauchy sequence, see Theorem 3.8. The convergence analysis is performed in Section 4.InTheorem4.1 we prove that the triple (un i,kn i,p n i)n has a unique limit, which is a solution of the variational formulation (2.3)–(2.4). We finally present some numerical tests in Section 5. These tests are realized with the software FreeFEM3D (see [12]) in meaningful situations, that agree with the expectations of our result. 2. Iterative scheme We shall at first describe the weak formulation of problem (1.1). We assume that αiand γiare bounded functions from the set of nonnegative real numbers R+onto R, and belong to W1,∞(Ωi), which satisfy ∀∈R+,δ 1≥αi()≥νand δ1≥γi()≥ν, (2.1) and ∀∈R+,|α i()|≤δ2and |γ i()|≤δ2,(2.2) where δ1,δ 2and νare positive constants. System (1.1) admits the following variational formulation: Find (ui,p i,k i)∈Xi×L2(Ωi)×W1,r(Ωi) such that, for all (vi,q i,ϕ i)∈Xi×L2(Ωi)×W1,r 0(Ωi), ai(ki;ui,vi)+bi(vi,p i)+κiΓ|ui−uj|(ui−uj)·vidτ=Ωi fi·vidx bi(ui,q i)=0,(2.3) and, ki=0onΓ i,k i=λ|ui−uj|2on Γ,and Ci(ki;ki,ϕ i)=Ωi αi(ki)|∇ui|2ϕidx,(2.4) where the forms ai(·;·,·),b i(·,·)andCi(·;·,·) are defined by ai(i;ui,vi)=Ωi α(i)∇ui:∇vidx, bi(vi,q i)=−Ωi qi∇·vidx, Ci(i;ki,ϕ )=Ωi γi(i)∇ki·∇ϕidx. Note that the bilinear forms aiand Ciin (2.3) depend on ki. Remark 2.1. Since ui∈Xithen its trace on Γ belongs to H 1 2 00(Γ). Thus by using definition of this space ui|Γ belongs to H1 2(Γ) and applying the Sobolev embedding from H1 2(Γ) into L3(Γ)d, we conclude that the integral Γ|ui−uj|(ui−uj)·vidτis well defined. AN ITERATIVE PROCEDURE TO SOLVE A COUPLED TWO-FLUIDS TURBULENCE MODEL 697 This formulation makes sense, as αi(ki)|∇ui|2∈W1,r(Ωi)whenui∈Xi. In Lewandowski [19]itisproved that this formulation admits at least a solution. We shall consider the following iterative procedure: Given (un i,p n i,kn i)∈Xi×L2(Ωi)×W1,r(Ωi),i=1,2, obtain (un+1 i,p n+1 i,kn+1 i)∈Xi×L2(Ωi)×W1,r(Ωi), such that ∀(vi,q i,ϕ i)∈Xi×L2(Ωi)×W1,r 0(Ωi), ai(kn i;un+1 i,∇vi)+bi(vi,p n+1 i)+κiΓ|un+1 i−un+1 j|(un+1 i−un+1 j)·vidτ=Ωi fi·vidτ, (2.5) and bi(un+1 i,q i)=0,(2.6) and kn+1 i=0onΓ i,(2.7) kn+1 i=λ|un+1 1−un+1 2|2on Γ,(2.8) Ci(kn i;kn+1 i,ϕ i)=Ωi αi(kn i)|∇un+1 i|2ϕidx.(2.9) 3. Contractiveness In this section we prove that the sequence of TKE (kn i)nis contracting, and that consequently the sequences of velocities (un i)nalso is contracting, in the sense that ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 2  i=1 ∇(un+1 i−un i)2 L2(Ωi)≤K 2  i=1 ∇(kn i−kn−1 i)2 L2(Ωi),and 2  i=1 ∇(kn+1 i−kn i)2 L2(Ωi)≤K 2  i=1 ∇(kn i−kn−1 i)2 L2(Ωi). (3.1) We may interpret these inequalities in the sense that the sequence of pairs (un i,kn i)nis contracting in the Hilbert space Xi×L2(Ωi). However, to simplify our derivation, we shall not explicitly use this space. Finally, we show that the pressures (pn i)nis a Cauchy sequence. We suppose from now on that the sequences (un i)nand (kn i)nverify the following hypothesis. Hypothesis 3.1. ∀n∈N,un i∈W1,3+ε(Ωi)dand kn i∈W1,3(Ωi), and one has un iW1,3+ε(Ωi)d≤M, kn iW1,3(Ωi)≤M, where Mand εare two fixed positive numbers. Remark 3.2. Note that the natural estimates for velocities in model (1.1)areinH1norm, not in W1,3norm. Indeed, choosing viequal to 1 κiun+1 i∈Xiin equation (2.5), and summing upon i=1,2gives 2  i=1 1 κiΩi αi(kn i)|∇un+1 i|2dx+Γ|un+1 1−un+1 2|3dτ= 2  i=1 1 κiΩi fiun+1 idx. Since the integrated term on Γ is nonnegative and thanks to (2.1), we deduce ν cM 2  i=1 ∇un+1 i2 L2(Ωi)d≤1 cm 2  i=1 Ωi fiun+1 idx. 698 T. CHAC ´ ON REBOLLO ET AL. Using the Cauchy-Schwarz and Poincar´e-Friedrichs inequalities, we obtain 2  i=1 ∇un+1 i2 L2(Ωi)d≤c ν2 2  i=1 fi2 L2(Ωi)d,(3.2) where cis a positive constant, depending only on the domains Ωiand the friction coefficients κi. We next prove that the contractiveness of the TKE implies that of the velocities. Lemma 3.3. Assume that Hypothesis 3.1 holds and that fi∈L2(Ωi)d,i=1,2. Then there exists a positive constant c, depending only on Ωi, such that 2  i=1 ∇(un+1 i−un i)2 L2(Ωi)d≤cδ2 2M2 ν2 2  i=1 ∇(kn i−kn−1 i)2 L2(Ωi)d.(3.3) Proof. Let us take vi=1 κi (un+1 i−un i)∈Xias a test function in (2.5)atiterationsnand n+ 1. Then, calculating the difference between both obtained equations, and summing on i=1,2, yields 2  i=1 1 κiΩi αi(kn)∇un+1 i−un i:∇un+1 i−un idx + 2  i=1 1 κiΩiαi(kn)−αikn−1 i∇un i:∇un+1 i−un idx +Γun+1 1−un+1 2un+1 1−un+1 2−|un 1−un 2|(un 1−un 2)·un+1 1−un+1 2−(un 1−un 2)dτ=0. The following inequality holds for all vectors a,b∈Rd, (|b|b−|a|a)·(b−a)≥0.(3.4) To prove it, consider the function J:Rd→ Rdefined by J(a)=2 3|a|3.Jis convex and differentiable. Thus, (∇J(a)−∇J(b)) ·(b−a)≥0,∀a,b∈Rd. Then, (3.4) follows as ∇J(a)=|a|a·b. We deduce that 2  i=1 1 κiΩi αi(kn i)|∇(un+1 i−un i)|2dx+ 2  i=1 1 κiΩi (αi(kn i)−αi(kn−1 i))∇un i·∇(un+1 i−un i)dx≤0.(3.5) AN ITERATIVE PROCEDURE TO SOLVE A COUPLED TWO-FLUIDS TURBULENCE MODEL 699 It comes from Hypothesis 3.1 that ∇un ibelongs to L3(Ωi)dand that ∇un iL3(Ωi)d≤M. Furthermore, according to the relation (2.2) and the canonical injection from H1(Ωi)toL6(Ωi)andH¨older inequality, we obtain ν 2  i=1 ∇(un+1 i−un i)2 L2(Ωi)d≤δ2 2  i=1 Ωi|kn i−kn−1 i||∇un i||∇(un+1 i−un i)|dx ≤δ2 2  i=1 kn i−kn−1 iL6(Ωi)∇un iL3(Ωi)d∇(un+1 i−un i)L2(Ωi)d, ≤δ2 2M2 2ν 2  i=1 kn i−kn−1 i2 L6(Ωi)+ν 2 2  i=1 ∇(un+1 i−un i)2 L2(Ωi)d. From this estimate we conclude relation (3.3).  We next prove the contractiveness of the sequence of TKE (kn i)n. Theorem 3.4. Assume that Hypothesis 3.1 holds and that fi∈L2(Ωi)d. Then there exists a positive constant c, depending only on Ωiand on the data κiand λ, such that for all n∈N∗, 2  i=1 ∇(kn+1 i−kn i)2L2(Ωi)≤c(δ2 1+1)δ2 2 ν3M 2  i=1 ∇(kn i−kn−1 i)2L2(Ωi).(3.6) Proof. The proof of this theorem is made in several steps. First step. Choice of the test function We first choose a particular test function ϕiin the equations (2.7)–(2.9). For that purpose, we need to introduce the special space H 1 2 00(Γ) (see [20], Chap. 1, Thm. 11.7 for instance). We also need to introduce the following operator. Let Ribe a continuous harmonic lifting operator from H 1 2 00(Γ) to H1(Ωi), defined as follows. For any ηin H 1 2 00(Γ), Riηbelongs to H1(Ωi), and satisfies ⎧ ⎨ ⎩ −ΔRiη=0inΩ i, Riη=ηon Γ,and Riη=0onΓ i. Moreover, one has ∀η∈H 1 2 00(Γ),RiηH1(Ωi)≤cRηH 1 2 00(Γ),(3.7) where cR>0 depends only on Ωi. According to Hypothesis 3.1,∀n∈N∗,k n i∈W1,3(Ωi), then its trace on Γ belongs to W2 3,3(Γ). Thus, by Sobolev’s injections, it belongs to H1 2(Γ). Furthermore, kn i=0onΓ i,thenkn i|Γbelongs to H 1 2 00(Γ). The idea consists in choosing the test function ϕiequal to (kn+1 i−kn i)−Ri(kn+1 i−kn i)inequation(2.9)at steps nand n+ 1. Then, we make the difference between both obtained equations, and sum upon i=1,2. We find ν 2  i=1 ∇(kn+1 i−kn i)2 0,Ωi≤ 7  j=1 Ij,(3.8) 700 T. CHAC ´ ON REBOLLO ET AL. where I1= 2  i=1Ωi αi(kn i)(|∇un+1 i|2−|∇un i|2)(kn+1 i−kn i)dx , I2= 2  i=1Ωi (αi(kn i)−αi(kn−1 i))|∇un i|2(kn+1 i−kn i)dx , I3= 2  i=1Ωi (γi(kn i)−γi(kn−1 i))∇kn i·∇(kn+1 i−kn i)dx , I4= 2  i=1Ωi (γi(kn i)−γi(kn−1 i))∇kn i·∇Ri(kn+1 i−kn i)dx , I5= 2  i=1Ωi αi(kn i)(|∇un+1 i|2−|∇un i|2)Ri(kn+1 i−kn i)dx , I6= 2  i=1Ωi (αi(kn i)−αi(kn−1 i))|∇un i|2Ri(kn+1 i−kn i)dx ,and I7= 2  i=1Ωi γi(kn i)∇(kn+1 i−kn i)·∇Ri(kn+1 i−kn i)dx . Second step. Estimates of Ij,1≤j≤7 Estimation of I1.We write |∇un+1 i|2−|∇un i|2=∇un+1 i−un i·∇un+1 i+un i, and use Hypothesis 3.1,and relation (2.1). Thanks to the Sobolev embedding of H1(Ωi)intoL6(Ωi)andfromH¨older and Poincar´e-Friedrichs inequalities, we obtain Ωi αi(kn i)(|∇un+1 i|2−|∇un i|2)(kn+1 i−kn i)dx ≤δ1∇(un+1 i−un i)L2(Ωi)dun+1 iW1,3(Ωi)d+un iW1,3(Ωi)dkn+1 i−kn iL6(Ω) ≤Mcδ1∇(un+1 i−un i)L2(Ωi)d∇(kn+1 i−kn i)L2(Ωi), where cis a positive constant, depending only on domains Ωi. To simplify the calculations, we introduce a positive number βwhich we shall fix later. According to Young’s inequality 1 βa2+βb2≥2ab, ∀a, b ∈R,and ∀β>0,(3.9) we obtain Ωi αi(kn i)(|∇un+1 i|2−|∇un i|2)(kn+1 i−kn i)dx≤ν β∇(kn+1 i−kn i)2 L2(Ωi)+βM2δ2 1c2 ν∇(un+1 i−un i)2 L2(Ωi)d. AN ITERATIVE PROCEDURE TO SOLVE A COUPLED TWO-FLUIDS TURBULENCE MODEL 701 Summing upon i=1,2, and due to the relation (3.3) from Lemma 3.3, there exists a positive constant c1, depending only on Ωi,αiand M, such that I1≤c1βδ2 1δ2 2M2 ν3 2  i=1 ∇(kn i−kn−1 i)2 L2(Ωi)+ν β 2  i=1 ∇(kn+1 i−kn i)2 L2(Ωi).(3.10) Estimation of I2and I3.Using the same arguments we used for estimation of I1, there exists two positive constants, depending only on Ωi,γ iand M, such that I2≤c2βδ2 2M2 ν 2  i=1 ∇(kn i−kn−1 i)2 L2(Ωi)+ν β 2  i=1 ∇(kn+1 i−kn i)2 L2(Ωi),(3.11) and I3≤c3βδ2 2M2 ν 2  i=1 ∇(kn i−kn−1 i)2 L2(Ωi)+ν β 2  i=1 ∇(kn+1 i−kn i)2 L2(Ωi).(3.12) Estimation of I4.We recall that I4= 2  i=1Γi (γi(kn i)−γi(kn−1 i))∇kn i·∇Ri(kn+1 i−kn i)dx . Let us apply the Mean Value Theorem to the function γi,userelation(2.2)andH¨older inequality. We find Γi (γi(kn i)−γi(kn−1 i))∇kn i·∇Ri(kn+1 i−kn i)dx≤δ2kn i−kn−1 iL6(Ωi)∇kn iL3(Ωi)∇Ri(kn+1 i−kn i)L2(Ωi). The continuity of the lifting operator (3.7)andHypothesis3.1 imply Γi (γi(kn i)−γi(kn−1 i))∇kn i·∇Ri(kn+1 i−kn i)dx≤cRδ2Mkn i−kn−1 iL6(Ωi)kn+1 i−kn iH 1 2 00(Γ). According to the continuity of the canonical injection from H1(Ωi)toL6(Ωi), the continuity of the trace operator from Xito H 1 2 00(Γ), and using Young’s inequality (3.9), there exists a positive constant c4>0, depending only on Ωi,γ i,and M, such that Γi (γi(kn i)−γi(kn−1 i))∇kn i·∇Ri(kn+1 i−kn i)dx≤c4βδ2 2M2 ν∇(kn i−kn−1 i)2 L2(Ωi)+ν β∇(kn+1 i−kn i)2 L2(Ωi). Summing on i=1,2, we find I4≤c4βδ2 2M2 ν 2  i=1 ∇(kn i−kn−1 i)2 L2(Ωi)+ν β 2  i=1 ∇(kn+1 i−kn i)2 L2(Ωi).(3.13) Estimation of I5.We have I5= 2  i=1Ωi αi(kn i)(|∇un+1 i|2−|∇un i|2)Ri(kn+1 i−kn i)dx . 708 T. CHAC ´ ON REBOLLO ET AL. embedding from W1,r 0(Ωi)toH1 0(Ω), we can write for m≥m0 ϕm i−ϕiH1 0(Ωi)≤η 3·(4.3) Thanks to the triangular inequality, we obtain Ωi γi(kn i)∇kn+1 i·∇ϕidx−Ωi γi(ki)∇ki·∇ϕidx≤Ωi γi(kn i)∇kn+1 i·∇ϕidx−Ωi γi(kn i)∇kn+1 i·∇ϕm idx +Ωi γi(kn i)∇kn+1 i·∇ϕm idx−Ωi γi(ki)∇(ki)·∇ϕm idx +Ωi γi(ki)∇(ki)·∇ϕm idx−Ωi γi(ki)∇(ki)·∇ϕidx . (4.4) H¨older inequality implies Ωi γi(kn i)∇kn+1 i·∇ϕidx−Ωi γi(ki)∇ki·∇ϕidx≤δ1∇kn+1 iL2(Ωi)∇(ϕm i−ϕi)L2(Ωi) +Ωi γi(kn i)∇kn+1 i·∇ϕm idx−Ωi γi(ki)∇(ki)·∇ϕm idx +δ1∇ki0∇(ϕm i−ϕi)L2(Ωi). Using relations (4.2)–(4.3), and the strong convergence of the sequence (kn i)nin H1(Ωi), we deduce that for any function ϕi∈W1,r 0(Ωi), we have lim n→∞ Ωi γi(kn i)∇(kn+1 i)·∇ϕidx=Ωi γi(ki)∇(ki)·∇ϕidx. To show the second equation of the relation (4.1), we write Ωi αi(kn i)|∇un+1 i|2−αi(ki)|∇ui|2ϕidx ≤Ωi αi(kn i)|∇un+1 i|2−|∇ui|2|ϕi|dx+Ωi|αi(kn i)−αi(ki)||∇ui|2|ϕi|dx, thanks to Hypothesis 3.1, the Sobolev embedding from H1(Ωi)toL6(Ωi)andH¨older inequality, we can write for all ϕi∈D(Ωi), Ωi αi(kn i)|∇un+1 i|2−αi(ki)|∇ui|2ϕidx ≤δ1ϕi∞∇un+1 iL2(Ωi)d∇uiL2(Ωi)d∇(un+1 i−ui)L2(Ωi)d+cM2δ2ϕi∞kn i−kiH1(Ωi), where cis a positive constant, depending only on the Ωi. Due to the strong convergence in H1of the sequences (un i)nand (kn i)nto uiand ki, Ωiαi(kn i)|∇un+1 i|2−αi(ki)|∇ui|2ϕidx≤η 3· AN ITERATIVE PROCEDURE TO SOLVE A COUPLED TWO-FLUIDS TURBULENCE MODEL 709 Let now ϕi∈W1,r 0(Ωi). The density of D(Ωi)inW1,r 0(Ωi), implies that there exists a sequence (ϕm i)min D(Ωi), such that ϕm i−ϕiL3(Ωi)≤c(Ωi)ϕm i−ϕiW1,r (Ωi)≤η 3, where c(Ωi) only depends on Ωi. Finally, using the triangular inequality, we obtain Ωiαi(kn i)|∇un+1 i|2−αi(ki)|∇ui|2ϕidx≤Ωiαi(kn i)|∇un+1 i|2−αi(ki)|∇ui|2ϕm idx +Ωi αi(kn i)|∇un+1 i|2(ϕm i−ϕi)dx +Ωi αi(ki)|∇ui|2(ϕm i−ϕi)dx . Thus, Ωiαi(kn i)|∇un+1 i|2−αi(ki)|∇ui|2ϕidx≤αi∞∇un+1 i2 L3(Ωi)dϕm i−ϕiL3(Ωi) +ε 3+αi∞∇ui2 L3(Ωi)dϕm i−ϕiL3(Ωi). We deduce that ∀ϕi∈W1,r 0(Ωi),lim n→∞ Ωi αi(kn i)|∇un+1 i|2ϕidx=Ωi αi(ki)|∇ui|2ϕidx. Third step. Boundary conditions of TKE on Γ In this step, we show that ki=λ|u1−u2|2on Γ. Consider the following triangular inequality λ|u1−u2|2−kiH1 2(Γ) ≤λ|un 1−un 2|2−kiH1 2(Γ) +λ|un 1−un 2|2−λ|u1−u2|2H1 2(Γ). Due to the strong convergence of (kn i)nto kiin H1(Ωi), and thanks to the continuity of the trace operator from H1(Ωi)toH1 2(Γ), the sequence kn iconverges strongly to kiin H1 2(Γ). Furthermore, kn i=λ|un+1 1−un+1 2|2 on Γ. Thus, lim n→∞ λ|un+1 1−un+1 2|2=ki, strongly in H1 2(Γ). Let us now prove that lim n→∞ λ|un 1−un 2|2−λ|u1−u2|2H1 2(Γ) =0. Using the identity a2−b2=(a−b)(a+b) and Lemma 3.5,weobtain λ|un 1−un 2|2−λ|u1−u2|2H1 2(Γ) ≤λ||(un 1−u1)−(un 2−u2)||H1 2(Γ) || (un 1+u1)−(un 2+u2)||W1−1 3+ε,3+ε(Γ)d, so, λ|un 1−un 2|2−λ|u1−u2|2H1 2(Γ) ≤λun 1−u1H1 2(Γ) +un 2−u2H1 2(Γ)un 1+u1W1−1 3+ε,3+ε(Γ)d+un 2+u2W1−1 3+ε,3+ε(Γ)d. Due to the continuity of the trace operators from W1,3+ε(Ω1)dto W1−1 3+ε,3+ε(Γ)dand from H1(Ωi)toH1 2(Γ), and by Hypothesis 3.1, there exists a positive constant c, that only depends on the domains Ωi, such that λ|un 1−un 2|2−λ|u1−u2|2H1 2(Γ) ≤cMλ un 1−u1H1(Ω1)+un 2−u2H1(Ω2). 710 T. CHAC ´ ON REBOLLO ET AL. Finally, due to the strong convergence of (un i)nto uiin H1(Ωi), we deduce that lim n→∞ λ|un 1−un 2|2−λ|u1−u2|2H1 2(Γ) =0. Consequently, ki=λ|u1−u2|2on Γ. This, finishes the proof of this theorem.  5. Numerical experiments To conclude this paper, we use the algorithm introduced to solve the interaction of the ocean and the atmosphere in a simplified geometry. The discretization is performed using a Spectral method based on Legendre polynomials (see [2,10]or[11] for instance) that we have implemented in FreeFEM3D4. The algorithm presented in this paper is nonlinear, thus it cannot be used “as is”. Various strategies could be employed to treat this nonlinearity (Newton or fixed point algorithms, etc.). According to the monotonic nature of the nonlinear friction boundary condition on the boundary Γ, we choose to linearize the term Γ|un+1 i−un+1 j|(un+1 i−un+1 j)·vidτ of Problem 1. We replace it by Γ|un i−un j|(un+1 i−un+1 j)·vidτ. This term is linear and still monotonic in the unknowns (un+1 1,un+1 2), due to property (3.4). The problems for (un+1 1,p n+1 1)and(un+1 2,p n+1 2) still are coupled but the overall problem is linear and admits a unique solution: Obtain (un+1 i,p n+1 i,kn+1 i)∈Xi×L2(Ωi)×W1,r(Ωi), such that ∀(vi,q i,ϕ i)∈Xi×L2(Ωi)×W1,r 0(Ωi), ai(kn i;un+1 i,∇vi)+bi(vi,p n+1 i)+κiΓ|un i−un j|(un+1 i−un+1 j)·vidτ=Ωi fi·vidτ, (5.1) and bi(un+1 i,q i)=0.(5.2) To perform a somewhat realistic computation, we consider turbulent viscosities αiand γiwith the structure νt+√k. We consider the data reported in [5]: •Geometry: –Ω1=]0,5[ ×]0,1[ ×]0,1[ describes the atmosphere; –Ω1=]0,5[ ×]0,−1[ ×]0,1[ is the ocean. •Physical data (taken from [5]): –γ1(k1)=3×10−3+0.277 ×10−4√k1; –γ2(k2)=3×10−2+0.185 ×10−5√k2; –αi(·)=γi(·), for i=1,2. •Friction coefficients (coming also from [5]): –κi=10 −3,fori=1,2; and –λ=5×10−2. These data correspond to an air-sea flow, each modeled by a simplified TKE-mixing layer turbulence model. The mixing lengths are calculated by wall laws. The friction coefficients, in their turn, are just tentative. The physical units are in MKS system. The velocity boundary conditions imposes that u1=0on Γ1\˜ Γ1, u1=(1,0,0) on ˜ Γ1and u2=0on Γ2,where˜ Γ1is the upper face (y=1)ofΩ 1. These settings are chosen in order to create a driven cavity-like flow in Ω1. One expects to generate another driven cavity-like flow in Ω2 rotating in the opposite sense. 4http://www.freefem.org/ff3d/. AN ITERATIVE PROCEDURE TO SOLVE A COUPLED TWO-FLUIDS TURBULENCE MODEL 711 (a) u1on the plane z=0.5(b)u1on the plane x=2.5 (c) u2on the plane z=0.5(d)u2on the plane x=2.5 Figure 1. Velocity fields on cutting planes. (a) k1on the plane z=0.5(b)k2on the plane z=0.5 Figure 2. Turbulent Kinetic Energy on cutting plane. The results are obtained using a (PN)3×PN−2discretization of (ui,p i) to avoid spurious modes of the Stokes problem. The TKE is discretized using a PNspace. For this particular simulation we have chosen the following degrees for uiand ki: 28 in direction xand 8 in directions yand z. Computed velocity fields uiand the TKE kiare represented in Figures 1and 2. The results are quantitatively correct: The atmosphere flow generates a driven-lid like flow in the ocean, due to the boundary conditions on the velocity at z= 0. Also, there is generation of TKE at z= 0, again due to the TKE generation boundary conditions. 712 T. CHAC ´ ON REBOLLO ET AL. 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1 0 1 2 3 4 5 6 7 8 9 u1 u2 k1 k2 Figure 3. Convergence history: computed L2norm of the difference of successive iterates. Figure 3shows the expected exponential convergence rate of the algorithm due to its contractiveness. Note that this does not depend on the type of discretization (see [24] where Finite Element approximations have also been used). 6. Conclusion In this paper, we have presented and analyzed a numerical scheme for the approximation of a model of two steady turbulent fluids with coupling at the interface. This is a simplified model for the atmosphere-ocean interaction, where we have neglected Coriolis forces and buoyancy effects, but have kept several non-linear interactions across the common boundary. The proposed scheme is mainly linear, a monotone nonlinearity being just kept at the interface between the fluids. We showed the convergence of the triple (un i,p n i,kn i), for reasonable hypothesis on the regularity of the velocity and the turbulent kinetic energy. This contribution ends with some numerical results, in good agreement with the theoretical expectations. Notably, the exponential convergence that appeared through the contractiveness of the sequences (un i)nand (kn i)nis found in our tests. Several extensions to this work can be considered. To name a few, taking into consideration anisotropic diffusion and Coriolis forces (more realistic) are straightforward generalizations of the present analysis. Also, switching to unsteady incompressible flows should be possible using our approach. Taking into account buoyancy effects is more technically involved, but it should also be possible, similarly to the extension of the standard analysis for incompressible Navier–Stokes to buoyancy effects. Acknowledgements. The authors are grateful to Christine Bernardi for her help and encouragements. References [1] J.J.F. Adams and R.A. Fournier, Sobolev spaces. Second edition, Pure and Applied Mathematics Series, Elsevier/Academic Press (2003). [2] C. Bernardi and Y. Maday, Approximations spectrales de probl`emes aux limites elliptiques,Mathematics & Applications 10. Springer-Verlag (1992). [3] C. Bernardi, T. Chac´on Rebollo, R. Lewandowski and F. Murat, A model for two coupled turbulent fluids. I. Analysis of the system, in Nonlinear partial differential equations and their applications, Coll`ege de France Seminar,Vol.XIV (Paris, 1997/1998), Stud. Math. Appl. 31, Amsterdam, North-Holland (2002) 69–102. AN ITERATIVE PROCEDURE TO SOLVE A COUPLED TWO-FLUIDS TURBULENCE MODEL 713 [4] C. Bernardi, T. Chac´on Rebollo, R. Lewandowski and F. Murat, A model for two coupled turbulent fluids. II. Numerical analysis of a spectral discretization. SIAM J. Numer. Anal. 40 (2003) 2368–2394. [5] C. Bernardi, T. Chac´onRebollo,M.G´omez M´armol, R. Lewandowski and F. Murat, A model for two coupled turbulent fluids. III. Numerical approximation by finite elements. Numer. Math. 98 (2004) 33–66. [6] C. Bernardi, T. Chac´on Rebollo, F. Hecht and R. Lewandowski, Automatic insertion of a turbulence model in the finite element discretization of the Navier-Stokes Equations. Math. Mod. Meth. Appl. Sci. 19 (2009) 1139–1183. [7] H. Brezis, Analyse Fonctionnelle : Th´eorie et Applications. Collection “Math´ematiques Appliqu´ees pour la Maˆıtrise”, Masson (1983). [8] F. Brossier and R. Lewandowski, Impact of the variations of the mixing length in a first order turbulent closure system. ESAIM: M2AN 36 (2002) 345–372. [9] K. Bryan, A numerical method for the study of the circulation of the world ocean. J. Comput. Phys. 4(1969) 347–369. [10] C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral methods – Fundamentals in single domains. Springer, Berlin, Germany (2006). [11] C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral methods – Evolution to complex geometries and applications to fluid dynamics. Springer, Berlin, Germany (2007). [12] S. Del Pino and O. Pironneau, A fictitious domain based on general pde’s solvers, in Proc. ECCOMAS 2001, Swansea, K. Morgan Ed., Wiley (2002). [13] V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations, Theory and Algorithms. Springer-Verlag, Germany (1986). [14] P. Grisvard, Elliptic Problems in Nonsmooth Domains,Monographs and Studies in Mathematics 24. Pitman (Advanced Publishing Program), Boston, USA (1985). [15] E. Hebey, Nonlinear analysis on manifolds: Sobolev spaces and inequalities,Courant Lecture Notes 5. American Mathematical Society, USA (1999). [16] B.E. Launder and D.B. Spalding, Mathematical Modeling of Turbulence. Academic Press, London, UK (1972). [17] J. Lederer and R. Lewandowski, A RANS 3D model with unbounded eddy viscosities. Ann. Inst. H. Poincar´eAnal.Non Lin´eaire 24 (2007) 413–441. [18] R. Lewandowski, Analyse Math´ematique et Oc´eanographie. Collection Recherches en Math´ematiques Appliqu´ees, Masson (1997). [19] R. Lewandowski, The mathematical analysis of the coupling of a turbulent kinetic energy equation to the Navier-Stokes equation with an eddy viscosity. Nonlinear Anal. 28 (1997) 393–417. [20] J.-L. Lions and E. Magenes, Probl`emes aux limites non homog`enes et applications 3,Travaux et Recherches Math´ematiques 20. Dunod, Paris, France (1970). [21] B. Mohammadi and O. Pironneau, Analysis of the k-epsilon turbulence model. RAM: Research in Applied Mathematics. Masson, Paris (1994). [22] J. Piquet, Turbulent Flows, Models and Physics. Springer, Germany (1999). [23] D.C. Wilcox, Turbulence Modeling for CFD. Sixth edition, DCW Industries, inc. California, USA (2006). [24] D. Yakoubi, Analyse et mise en œuvre de nouveaux algorithmes en m´ethodes spectrales.Ph.D.Thesis,Universit´e Pierre et Marie Curie, Paris, France (2007).