Time-averaged shallow water model: Asymptotic derivation and numerical validation

J. Math. Anal. Appl. 428 (2015) 930–950

Contents lists available at ScienceDirect

Journal of Mathematical Analysis and Applications www.elsevier.com/locate/jmaa

Time-averaged shallow water model: Asymptotic derivation and numerical validation ✩ J.M. Rodríguez a,∗ , R. Taboada-Vázquez b a

Department of Mathematics and Representation Methods, College of Architecture, University of A Coruña, Campus da Zapateira, 15071, A Coruña, Spain b Department of Mathematics and Representation Methods, School of Civil Engineering, University of A Coruña, Campus de Elviña, 15071, A Coruña, Spain

a r t i c l e

i n f o

Article history: Received 28 May 2014 Available online 20 March 2015 Submitted by W. Layton Keywords: Shallow waters Asymptotic analysis Reynolds Averaged Navier–Stokes equations (RANS) Large Eddy Simulation (LES) Modeling Filtering

a b s t r a c t The objective of this paper is to derive, from the Navier–Stokes equations in a shallow domain, a new bidimensional shallow water model able to filter the high frequency oscillations that are produced, when the Reynolds number is increased, in turbulent flows. With this aim, the non-dimensional Navier–Stokes equations are time-averaged, and then asymptotic analysis techniques have been used as in our previous works (Rodríguez and Taboada-Vázquez, 2005–2012 [8–14]). The small non-dimensional parameter considered, ε, is the quotient between the typical depth of the basin and the typical horizontal length of the domain; and it is studied what happens when ε becomes small. Once the new model has been justified, by the method of asymptotic expansions, we perform some numerical experiments. The results of these experiments confirm that this new model is able to approximate analytical solutions of Navier–Stokes equations with more accuracy than classical shallow water models, when high frequency oscillations appear. To reach a given accuracy, the time step for the new model can be much larger (even four hundred times larger) than the time step required for the classical models. © 2015 Elsevier Inc. All rights reserved.

1. Introduction As it is well known, the equations governing the behavior of a fluid are the Navier–Stokes equations. Due to their strong nonlinearity, high frequency oscillations are produced when the Reynolds number is increased, and the flow becomes unstable and turbulent. It is computationally very expensive to solve the equations directly, so at the moment, the most common approach in hydraulic engineering practice is to ✩ This work has been partially supported by Ministerio de Economía y Competitividad of Spain under grant MTM2012-36452C02-01 with the participation of FEDER. * Corresponding author. E-mail addresses: [email protected] (J.M. Rodríguez), [email protected] (R. Taboada-Vázquez).

http://dx.doi.org/10.1016/j.jmaa.2015.03.050 0022-247X/© 2015 Elsevier Inc. All rights reserved.

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

931

solve the Reynolds Averaged Navier–Stokes equations, in which the effect of turbulence is modeled rather than solved. We approximate Navier–Stokes equations using a shallow water model, but if the flow is turbulent, a very small time step must be chosen. This paper is focused on the derivation of a new bidimensional time-averaged shallow water model able to reduce these oscillations, and then able to achieve good results with larger time steps. Filtering has given good results when working with turbulent Navier–Stokes equations (see [15]). In the literature, we can find that the separation between large and small scales is traditionally assumed to be obtained by applying a spatial filter to the Navier–Stokes equations (see [1,4]); but time filtering is also suggested by several authors (see [2,7]). In this work, we shall use a time filter, thus avoiding to model the spacial filtered stress tensor (see [1]). Asymptotic analysis has been applied successfully to derive and justify shallow water models. The new model, developed from the incompressible Navier–Stokes equations with free surface, has been deduced in the spirit of the method proposed in our previous works [8–14] and [16]. In order to obtain a shallow water model, we consider a domain with small depth compared with its other dimensions. We use in the sequel the thin-layer assumption and introduce a “small” non-dimensional parameter ε = HC /LC where HC and LC are, respectively, the typical scales for the vertical and the horizontal dimensions of the fluid domain of interest. The outline of the paper is as follows. In the next section we introduce and render non-dimensional the model that serves as our starting point. Then, the time averaging process is described in Section 3 and asymptotic analysis is applied following the ideas of [8–14] (Section 4) to derive our shallow water model that is presented in Section 5. We show, in Section 6, that this new model is able to obtain a given accuracy using time steps larger than the time steps needed by classical shallow water models. Finally, we make some concluding remarks in Section 7. 2. The three-dimensional model equations In this section we present the three-dimensional incompressible Navier–Stokes model that serves as the starting point for our subsequent development. The first subsection gives the basic mass and momentum balance laws for a basin with varying bottom topography and a free top surface, and supplements them with appropriate boundary conditions. In Section 2.2 we introduce the shallow water scaling, define a non-dimensional parameter and non-dimensionalize the three-dimensional model in terms of that parameter. 2.1. Three-dimensional incompressible flow Let us start with the Navier–Stokes system [5] for incompressible homogeneous fluids, with gravity and Coriolis force, evolving in a sub-domain of R3 . As the domain, the functions and variables involved in this problem depend on ε, we indicate this dependence with superscript ε. Therefore, we have the following general formulation expression:  ε  ∂U  ε · ∇ε U ε  ε = − 1 ∇ε P ε + νΔε U ε+F + U e ∂tε ρ0

(1)

ε=0 div U

(2)

and we consider this system for tε ∈ [0, T ], where:

(xε , y ε ) ∈ D ⊂ R2 ,

B ε (xε , y ε ) ≤ z ε ≤ S ε (tε , xε , y ε )

932

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

Fig. 1. Notations: water height H ε (tε , xε , y ε ), free surface S ε (tε , xε , y ε ) and bottom B ε (xε , y ε ).

• S ε (tε , xε , y ε ) represents the free surface elevation (unknown) and B ε the bathymetry (it is not constant and it is supposed to be known). The water height is H ε = S ε − B ε (see Fig. 1)  ε = (U ε (tε , xε , y ε , z ε ), U ε (tε , xε , y ε , z ε ), U ε (tε , xε , y ε , z ε )) is the three-dimensional velocity of the fluid • U 1 2 3 • P ε (tε , xε , y ε , z ε ) is the pressure • ρ0 denotes the density of the fluid • ν is the kinematic viscosity  U  ε = −gk−2Φ×  ε is the volume force per unit mass, where g is the gravitational acceleration (assumed • F e  ×U  ε is the Coriolis acceleration (where the angular velocity of rotation of the Earth constant)and −2Φ   = Φ sin ϕk + cos ϕj with Φ = 7.29 × 10−5 rad/s; ı, j and k denote the unit vectors pointing East, is Φ North and vertically upward (respectively); ϕ, the North latitude, is considered constant). The kinematic continuity condition U3ε =

∂H ε ∂S ε ∂S ε + U1ε ε + U2ε ε ε ∂t ∂x ∂y

at z ε = S ε (tε , xε , y ε )

(3)

is rewritten in equivalent form (for incompressible fluids): ∂ ∂H ε + ∂tε ∂xε

S

ε

∂ U1ε dz ε + ε ∂y

Bε

S

ε

U2ε dz ε = 0

(4)

Bε

Eqs. (1)–(2) must be supplemented by boundary conditions. • At the bottom, – the non-penetration condition is satisfied:  ε · nε = 0 at z ε = B ε (xε , y ε ) U

(5)

where nε denotes the outward unit normal to the boundary of the domain; – tangential forces must be equal to the friction force: ε (I − nε ⊗ nε )Tε nε = F R

at z ε = B ε (xε , y ε )

(6)

 ε |U  ε (C ε is small),  ε = −ρ0 C ε |U Typically, the friction force per unit of surface area is of the form F R R R see for example [3] or [17].

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

933

• On the free top surface we assume that the only external force acting on the fluid is the wind stress. In particular, we assume that the surface tension and ambient atmospheric pressure variations are negligible. This leads to the boundary condition  ε at z ε = S ε (tε , xε , y ε ) Tε nε = −Psε (tε , xε , y ε )nε + F W

(7)

 ε is the force of the wind and P ε is the atmospheric pressure at the surface (supposed known). where F s W The stress tensor (Tε ) is given by:  Tijε = −P ε δij + μ

∂Ujε ∂Uiε + ε ∂xj ∂xεi

 i, j = 1, 2, 3

(8)

where xε1 = xε , xε2 = y ε , xε3 = z ε , μ = ρ0 ν is the dynamic viscosity and δij is Kronecker’s delta. We also suppose that the incoming and outcoming flows are known at each instant (other kind of boundary conditions may be considered). Finally, initial conditions must be imposed too. 2.2. Non-dimensionalization The shallow water approximation is characterized by the smallness of a non-dimensional parameter that we can identify assuming that the typical depth of the basin (HC ) is much smaller than the typical horizontal length (LC ), i.e., that HC =ε LC

where

ε << 1

(9)

This small parameter is an aspect ratio. We introduce the non-dimensional independent variables t, x, y and z by t=

tε , TC

x=

xε , LC

y=

yε , LC

z=

zε εLC

(10)

where TC is a typical time. Recalling (9), the non-dimensional water height and bottom surface are defined by h(t, x, y) =

H ε (tε , xε , y ε ) , εLC

b(x, y) =

B ε (xε , y ε ) εLC

(11)

so the non-dimensional top surface is given by s(t, x, y) =

S ε (tε , xε , y ε ) εLC

(12)

We shall now introduce non-dimensional functions and constants: uεi (t, x, y, z) =

TC ε ε ε ε ε U (t , x , y , z ), LC i

pε (t, x, y, z) =

TC2 P ε (tε , xε , y ε , z ε ), ρ0 L2C

G=

TC2 g, LC

φ = TC Φ

i = 1, 2, 3 ps (t, x, y) =

(13) TC2 P ε (tε , xε , y ε ) ρ0 L2C s

(14) (15)

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

934

ε σij (t, x, y, z) =

TC2 T ε (tε , xε , y ε , z ε ) ρ0 L2C ij

TC2  ε ε fW = F , ρ0 L2C W

i, j = 1, 2, 3

(16)

TC2  ε fRε = F ρ0 L2C R

(17)

where we have assumed that ps does not depend on ε. 2.2.1. The non-dimensional equations We now express the three-dimensional incompressible Navier–Stokes system (1), (2) and (4) in terms of the above non-dimensional variables and functions:   1 ∂uε ∂pε 1 ∂ 2 uε1 ∂ 2 uε1 1 ∂ 2 uε1 ∂uε1 ∂uε ∂uε + uε1 1 + uε2 1 + uε3 1 = − + + + ∂t ∂x ∂y ε ∂z ∂x Re ∂x2 ∂y 2 ε2 ∂z 2 + 2φ ((sin ϕ) uε2 − (cos ϕ) uε3 ) ∂uε2 ∂t

+

∂uε uε1 2 ∂x

+

∂uε uε2 ε2 ∂y

+

1 ε ∂uε2 u3 ε

∂z

=−

ε

∂p 1 + ∂y Re



1 ∂uε 1 ∂pε 1 ∂uε ∂uε ∂uε3 + uε1 3 + uε2 3 + uε3 3 = − + ∂t ∂x ∂y ε ∂z ε ∂z Re

∂ 2 uε2 ∂x2



+

∂ 2 uε2 ∂y 2

+

1 ∂ 2 uε2 ε2 ∂z 2

(18)



∂ 2 uε3 ∂ 2 uε3 1 ∂ 2 uε3 + + 2 2 2 ∂x ∂y ε ∂z 2

− 2φ (sin ϕ) uε1 

− G + 2φ (cos ϕ) uε1 ∂uε1

∂uε2

(20)

∂uε3

1 + =0 ∂y ε ∂z s s ∂ ∂ ∂h ε + u1 dz + uε2 dz = 0 ∂t ∂x ∂y ∂x

+

b

where ν

(19)

(21)

(22)

b

TC 1 . The non-penetration condition (5) yields: = L2C Re uε3

  ε ∂b ε ∂b + u2 at z = b = ε u1 ∂x ∂y

The non-dimensional stress tensor can be written:   ∂uεj 1 ∂uεi ε ε i, j = 1, 2 σij = −p δij + + Re ∂xj ∂xi   ∂uε3 1 1 ∂uεi ε σi3 i = 1, 2 + = Re ε ∂z ∂xi ε σ33 = −pε +

2 1 ∂uε3 Re ε ∂z

so boundary conditions (6)–(7) result ⎛ ∂b ⎞ ∂b ε ε ε ε σ11 + ε σ12 − σ13 ∂y ⎜ ∂x ⎟   ⎜ ⎟ 2 ε ⎜ ∂b ε ⎟ ∂b 2 1 ∂b ∂u1 ε ε ε ⎟ ⎜ ε σ12 + ε σ22 − σ23 ⎟ + p − 2  2 ε  ⎜ ∂x ∂y Re ∂x ∂x ∂b ∂b ⎜ ⎟ 1+ ε + ε ⎝ ∂b ⎠ ∂b ∂x ∂y ε ε ε ε σ13 + ε σ23 − σ33 ∂x ∂y  ε     2 ε ∂uε2 ∂b 1 ∂uε1 ∂uε3 ∂b 1 ∂uε3 ∂u2 ∂u1 2 ∂b ∂b + −ε + + ε + +ε ∂x ∂y ∂y ∂x ∂x ε ∂z ∂x ∂y ∂y ε ∂z

(23)

(24) (25) (26)

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

935

⎛ ∂b ⎞  ε ⎞ ⎛     ⎜ ∂x ⎟   2  fRε 1 2 ε ε  ⎜ ⎟ ∂u3 ∂b 1 ∂u2 ∂b ∂b ⎜ ∂b ⎟ = 1 + ε2 ⎝ fRε ⎠ + −ε + ⎜ε ⎟ 2 ∂y ε ∂z ∂y ∂x ∂y ⎝ ∂y ⎠ fRε 3 −1 at z = b

(27)

⎛

⎞

∂s ε ∂s ε ⎛ ε ∂s ⎞ −ε σ11 − ε σ12 + σ13 −ε ∂x ∂y ⎜ ⎟ ⎛ ε ⎞ ⎜ ⎟ ⎜ 2  2 fW  ∂x ⎟  ⎜ ⎜ ⎟ ⎟ ∂s ε ∂s ⎝ ε 1 ⎠ ε ε ⎟ ⎜ −ε ∂s σ12 ∂s ⎟ + 1 + ε ∂s − ε σ22 + σ23 = −ps ⎜ + ε fW2 ⎜ ⎟ ⎜ ⎟ −ε ∂x ∂y ∂x ∂y ⎜ ⎟ ε ⎝ ⎠ ∂y fW ⎝ ⎠ 3 ∂s ε ∂s ε ε 1 −ε σ13 − ε σ23 + σ33 ∂x ∂y at z = s

(28)

3. Time averaging process The formal derivation of our shallow water model has two steps. We first obtain the time-averaged Navier–Stokes equations and then make the shallow water approximation. We define: u ¯εi (t, x, y, z; η)

1 = 2η

t+η  uεi (r, x, y, z)dr

i = 1, 2, 3

(29)

t−η

and, in the same way, we also define the other averaged functions of the problem. If we assume that 0 < η << 1, and we approximate uεi by its Taylor’s expansion for r ∈ (t − η, t + η): uεi (r, x, y, z) = uεi (t, x, y, z) + +

1 ∂ 2 uεi ∂uεi (t, x, y, z)(r − t) + (t, x, y, z)(r − t)2 ∂t 2 ∂t2

1 ∂ 3 uεi (t, x, y, z)(r − t)3 + · · · 3! ∂t3

(30)

then we have: u ¯εi (t, x, y, z; η) = uεi (t, x, y, z) +

∂ 4 uεi ∂ 2 uεi η2 η4 (t, x, y, z) + (t, x, y, z) + · · · 2 4 ∂t 3! ∂t 5!

(31)

It is easy to demonstrate for regular enough functions f ε and g ε , and small enough η, that the following equalities are fulfilled: ∂f ε ∂ f¯ε (t, x, y, z; η) = (t, x, y, z; η) ∂t ∂t ∂ 2 f¯ε η2 f ε (t, x, y, z) = f¯ε (t, x, y, z; η) − (t, x, y, z; η) + O(η 4 ) ∂t2 6 f ε g ε (t, x, y, z; η) = f¯ε (t, x, y, z; η)¯ g ε (t, x, y, z; η) ∂ f¯ε ∂¯ gε η2 + (t, x, y, z, η) (t, x, y, z, η) + O(η 4 ) ∂t ∂t 3  ¯ε ε  2 2 ε ∂ g η ∂ f ∂g f ε g ε = f¯ε g ε + 2 + f¯ε 2 + O(η 4 ) 6 ∂t ∂t ∂t  f ε (t, x, y, s(t, x, y)) = f¯ε z=s at z = s

(32) (33)

(34) (35) (36)

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

936

Applying (29), (32) and (34) to Eqs. (18)–(21), they yield:  ε 2 ε 2  ε 2 ε 2 ∂u ¯ε1 ¯ε1 ¯1 η ¯ε1 ¯1 η ∂u ¯1 ∂ u ∂u ¯2 ∂ u ε ∂u ε ∂u +u ¯1 + +u ¯2 + ∂t ∂x ∂t ∂t∂x 3 ∂y ∂t ∂t∂y 3   ε 2 ε 2 ε ∂u ¯3 ∂ u 1 ε ∂u ¯ ¯1 η ∂ p¯ε u ¯3 1 + + O(η 4 ) = − + ε ∂z ∂t ∂t∂z 3 ∂x   2 ε ∂2u 1 ∂2u 1 ∂ u ¯1 ¯ε1 ¯ε1 + 2φ ((sin ϕ) u ¯ε2 − (cos ϕ) u + + 2 ¯ε3 ) + 2 2 Re ∂x ∂y ε ∂z 2  ε 2 ε 2  ε 2 ε 2 ∂u ¯ε ¯2 η ∂u ¯ε ¯2 η ∂u ¯ε2 ∂u ¯1 ∂ u ∂u ¯2 ∂ u +u ¯ε1 2 + +u ¯ε2 2 + ∂t ∂x ∂t ∂t∂x 3 ∂y ∂t ∂t∂y 3   ε 2 ε 2 ∂u ¯3 ∂ u ∂ p¯ε ¯ε2 ¯2 η 1 ε ∂u u ¯3 + + O(η 4 ) = − + ε ∂z ∂t ∂t∂z 3 ∂y   2 ε 2 ε 2 ε ∂ u 1 ∂ u 1 ∂ u ¯2 ¯2 ¯2 − 2φ (sin ϕ) u ¯ε1 + + 2 + Re ∂x2 ∂y 2 ε ∂z 2  ε 2 ε 2  ε 2 ε 2 ∂u ¯1 ∂ u ∂u ¯2 ∂ u ¯ε3 ¯3 η ¯ε3 ¯3 η ∂u ¯ε3 ε ∂u ε ∂u +u ¯1 + +u ¯2 + ∂t ∂x ∂t ∂t∂x 3 ∂y ∂t ∂t∂y 3   ε 2 ε 2 1 ε ∂u 1 ∂ p¯ε ¯ε ¯3 η ∂u ¯3 ∂ u u ¯3 3 + + O(η 4 ) = − −G + ε ∂z ∂t ∂t∂z 3 ε ∂z   ∂2u 1 ∂2u 1 ∂2u ¯ε3 ¯ε3 ¯ε3 + 2φ (cos ϕ) u ¯ε1 + + + Re ∂x2 ∂y 2 ε2 ∂z 2 ∂u ¯ε2 1 ∂u ¯ε3 ∂u ¯ε1 + + =0 ∂x ∂y ε ∂z

(37)

(38)

(39) (40)

To average Eq. (22), in first place we write uε1 and uε2 using (33) and we get ∂ ∂h + ∂t ∂x

  s  s  ∂ ∂2u ¯ε1 η 2 ∂2u ¯ε2 η 2 ε u ¯ε1 − u ¯ dz + dz + O(η 4 ) = 0 − 2 ∂t2 6 ∂y ∂t2 6 b

(41)

b

This averaging process provides from (23)–(27):  ∂b ε ∂b +u ¯2 at =ε ∂x ∂y  ε 1 ∂u ¯i ε ε σ ¯ij = −¯ p δij + + Re ∂xj   ¯εi ∂u ¯ε3 1 1 ∂u ε σ ¯i3 + = Re ε ∂z ∂xi 

u ¯ε3

u ¯ε1

z=b ∂u ¯εj ∂xi

(42)

 i, j = 1, 2

(43)

i = 1, 2

(44)

¯ε3 2 1 ∂u ε σ ¯33 = −¯ pε + Re ε ∂z ⎛ ∂b ε ⎞ ∂b ε ε ε ∂x σ ¯11 + ε ∂y σ ¯12 − σ ¯13   2 ε ⎜ ∂b ε ⎟ 2 1 ∂b ∂u ¯1 ε ∂b ε ε ⎟ ⎜ε σ ¯22 − σ ¯23 ⎠ + p¯ − 2  2 ε  ⎝ ∂x ¯12 + ε ∂y σ Re ∂x ∂x ∂b ∂b ∂b ε ∂b ε ε 1+ ε + ε ε ∂x σ ¯13 + ε ∂y σ ¯23 − σ ¯33 ∂x ∂y ∂b ∂b +ε ∂x ∂y 2



∂u ¯ε2 ∂u ¯ε1 + ∂y ∂x



∂b −ε ∂x



¯ε1 ∂u ¯ε3 1 ∂u + ε ∂z ∂x





∂b + ε ∂y

2

∂u ¯ε2 ∂y

(45)

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

⎛ ∂b ⎞ ⎛ ⎞  ⎜ ε ∂x ⎟     2  2 f¯Rε 1 ε ε ε ⎟ ⎜ ¯2 ∂u ¯3 1 ∂u ¯3 ∂b ∂b 1 ∂ u ⎜ ¯ε ⎟ ⎜ ∂b ⎟ = 1 + ε2 ∂b + + −ε + ε2 ⎝ fR2 ⎠ ⎟ ⎜ε ∂y ε ∂z ∂y ε ∂z ∂x ∂y ⎝ ∂y ⎠ f¯Rε 3 −1

937

at z = b

(46) Condition (28) is time-averaged too and gives and similar but more complex expression, because s depends on t and we need to use formula (35). 4. Asymptotic derivation of the shallow water model After obtaining the time-averaged Navier–Stokes equations, we make the shallow water approximation by formally expanding solutions of the model in powers of ε. 4.1. Asymptotic expansion Our shallow water model derive from the assumption that the aspect ratio ε is small. We make the shallow water approximation by formally expanding solutions of the model in powers of ε. We seek a formal solution in the form of an asymptotic series: ⎧ ¯0i + ε¯ u1i + ε2 u ¯2i + · · · i = 1, 2, 3 u ¯ε = u ⎪ ⎪ ⎪ i ⎪ ⎪ ε 0 ⎪ p1 + ε2 p¯2 + · · · ⎪ p¯ = p¯ + ε¯ ⎨ ε 0 1 2 σ ¯ij =σ ¯ij + ε¯ σij + ε2 σ ¯ij + · · · i, j = 1, 2, 3 ⎪ ⎪ ⎪ ⎪ f¯Rε i = εf¯R1 i + ε2 f¯R2 i + · · · i = 1, 2 ⎪ ⎪ ⎪ ⎩ ¯ε 1 2 + ε2 f¯W + · · · i = 1, 2 fWi = εf¯W i i

(47)

ε The developments for f¯Rε i , f¯W (i = 1, 2) may begin in the term of order 1 because of their small order of i magnitude (see [16] for a rigorous justification). We substitute now this expansion into the system of Eqs. (37)–(46), and grouping the terms multiplied by the same power of ε, we arrive at a series of equations that will allow us to determine u ¯0i , p0 , etc. Let us show some particular examples. The value of u ¯03 can be found from the incompressibility condition (40) written at the leading order O(ε−1 ):

∂u ¯03 =0⇒u ¯03 = u ¯03 (t, x, y) ∂z

(48)

We now consider the boundary condition (42); to the leading order it becomes: u ¯03 = 0 for z = b

(49)

u ¯03 = 0

(50)

so

In the same way, upon inserting the above expansion into (44), to order O(ε−1 ), one finds that ∂u ¯0i = 0 i = 1, 2 ∂z

(51)

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

938

Now, incompressibility equation (40), to the leading order O(1), gives z  u ¯13

=−



∂u ¯01 ∂u ¯02 + ∂x ∂y

dz  + u ¯13 (t, x, y, b)

(52)

b

and the boundary condition (42), together with the fact that u ¯01 and u ¯02 do not depend on z ((51)), allows us to write:  u ¯13

= (b − z)

∂u ¯02 ∂u ¯01 + ∂x ∂y

 +u ¯01

∂b ∂b +u ¯02 ∂x ∂y

(53)

Eq. (39), to order O(ε−1 ), together with above expression for u ¯13 , leads us to conclude that p¯0 is constant with respect to the vertical variable: 1 ∂2u ∂ p¯0 ¯13 = =0 ∂z Re ∂z 2

(54)

and with the information that we obtain from boundary conditions we have: 2 p¯ = p¯s − Re 0



∂u ¯01 ∂u ¯02 + ∂x ∂y

 (55)

Eq. (41), to the leading order O(1), together with (51) becomes       ∂ ∂ η2 ∂ 2 u η2 ∂ 2 u ∂h ¯01 ¯02 0 + = O(η 4 ) + h u ¯01 − h u ¯ − 2 ∂t ∂x 6 ∂t2 ∂y 6 ∂t2

(56)

In summary, we are able to derive the following first terms and equations: u ¯03 = 0



∂u ¯k−1 ∂u ¯k−1 1 2 + ∂x ∂y  0  ∂u ¯02 2 ∂u ¯1 + p¯0 = p¯s − Re ∂x ∂y u ¯k3 = (b − z)

p¯ = (s − z)(G − 1

2φ(cos ϕ)¯ u01 )

(57)

 +u ¯k−1 1

∂b ∂b +u ¯k−1 , 2 ∂x ∂y

k = 1, 2

(58) (59)

2 − Re



∂u ¯12 ∂u ¯11 + ∂x ∂y



0 0 =σ ¯23 =0 σ ¯13

∂u ¯k2 ∂u ¯k1 = = 0 k = 0, 1 ∂z ∂z ∂u ¯01 η2 ∂ u ¯01 ∂ 2 u η2 ∂ u ¯02 ∂ 2 u ∂u ¯0 ¯01 ∂u ¯0 ¯01 +u ¯01 1 + +u ¯02 1 + + O(η 4 ) ∂t ∂x 3 ∂t ∂t∂x ∂y 3 ∂t ∂t∂y   1 ∂2u ∂2u ∂2u ¯01 ¯01 ¯21 ∂ p¯0 + + 2φ (sin ϕ) u ¯02 + + =− ∂x Re ∂x2 ∂y 2 ∂z 2  0 2 1  ∂u ¯11 η2 ∂ u ¯1 ∂ u ∂u ¯11 ∂ 2 u ∂u ¯1 ∂u ¯0 ¯1 ¯01 ∂u ¯1 ∂u ¯0 +u ¯01 1 + u ¯11 1 + + +u ¯02 1 + u ¯12 1 ∂t ∂x ∂x 3 ∂t ∂t∂x ∂t ∂t∂x ∂y ∂y  0 2 1  2 1 2 0 1 2 1 2 1 η ∂u ¯2 ∂ u η ∂u ¯3 ∂ u ¯1 ¯1 ∂u ¯ ¯1 ∂u ¯2 ∂ u + + +u ¯13 1 + + O(η 4 ) 3 ∂t ∂t∂y ∂t ∂t∂y ∂z 3 ∂t ∂t∂z

(60) (61) (62)

(63)

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

1 ∂ p¯1 + =− ∂x Re



∂2u ¯11 ∂2u ¯11 ∂2u ¯31 + + ∂x2 ∂y 2 ∂z 2



  + 2φ (sin ϕ) u ¯12 − (cos ϕ) u ¯13

∂u ¯02 ∂u ¯0 η2 ∂ u ¯01 ∂ 2 u ¯02 ∂u ¯0 η2 ∂ u ¯02 ∂ 2 u ¯02 +u ¯01 2 + +u ¯02 2 + + O(η 4 ) ∂t ∂x 3 ∂t ∂t∂x ∂y 3 ∂t ∂t∂y   ∂ p¯0 1 ∂2u ∂2u ∂2u ¯02 ¯02 ¯22 − 2φ (sin ϕ) u ¯01 =− + + + ∂y Re ∂x2 ∂y 2 ∂z 2  0 2 1  ∂u ¯12 ¯12 ¯02 ¯2 ¯02 ∂u ¯1 ∂u ¯0 η2 ∂ u ¯1 ∂ u ∂u ¯11 ∂ 2 u 0 ∂u 1 ∂u +u ¯1 +u ¯1 + + +u ¯02 2 + u ¯12 2 ∂t ∂x ∂x 3 ∂t ∂t∂x ∂t ∂t∂x ∂y ∂y  0 2 1  2 1 2 0 1 2 1 2 1 ∂u ¯2 ∂ u ¯2 ∂u ¯2 ∂ u ¯2 ∂u ¯ η ∂u ¯3 ∂ u ¯2 η + + +u ¯13 2 + + O(η 4 ) 3 ∂t ∂t∂y ∂t ∂t∂y ∂z 3 ∂t ∂t∂z   1 ∂2u ¯12 ∂2u ¯12 ∂2u ¯32 ∂ p¯1 + + + − 2φ (sin ϕ) u ¯11 =− ∂y Re ∂x2 ∂y 2 ∂z 2   2 0  2 0  u01 ) ∂(h¯ u02 ) η 2 ∂ ∂ u ∂ u ∂h ∂(h¯ ¯ ¯ ∂ + + − h 21 + h 22 = O(η 4 ) ∂t ∂x ∂y 6 ∂x ∂t ∂y ∂t   2 1  2 1  ∂(h¯ u11 ) ∂(h¯ u12 ) η 2 ∂ ¯1 ∂ ¯ ∂ u ∂ u = O(η 4 ) + − h 2 + h 22 ∂x ∂y 6 ∂x ∂t ∂y ∂t   ∂u ¯kj 1 ∂u ¯ki k k σ ¯ij = −¯ p δij + + i, j = 1, 2, k = 0, 1, · · · Re ∂xj ∂xi   ∂u ¯k3 1 ∂u ¯k+1 k i σ ¯i3 = + i = 1, 2, k = 1, 2, · · · Re ∂z ∂xi 2 1 ∂u ¯k+1 k 3 σ ¯33 = −¯ pk + k = 0, 1, · · · Re ε ∂z    0  2 ∂u ¯0 ∂u ¯02 ∂b 1 ∂u ¯1 ∂u ¯02 ∂b 1 2 1+ + + − f¯R1 1 at z = b σ ¯13 = Re ∂x ∂y ∂x Re ∂y ∂x ∂y  0   0  1 ∂u ∂u ¯02 ∂b 2 ∂u ∂u ¯02 ∂b ¯1 ¯1 1 σ ¯23 = + + +2 − f¯R1 2 at z = b Re ∂y ∂x ∂x Re ∂x ∂y ∂y    0  2 ∂u ¯01 ∂u ¯02 ∂s 1 ∂u ¯1 ∂u ¯02 ∂s 1 2 + + + σ ¯13 = Re ∂x ∂y ∂x Re ∂y ∂x ∂y      3 η2 ∂ ∂u ¯0 ∂u ¯02 ∂ 2 h ∂u ¯0 ∂u ¯02 ∂ h + 4 2 1+ +2 2 1 + 6Re ∂t ∂x ∂y ∂t∂x ∂x ∂y ∂t2 ∂x  0   0  3  ∂ ∂u ¯1 ∂u ¯02 ∂ 2 h ∂u ¯1 ∂u ¯02 ∂ h +2 + + + ∂t ∂y ∂x ∂t∂y ∂y ∂x ∂t2 ∂y

1 σ ¯23

1 + O(η 4 ) + f¯W at z = s 1  0   0  ¯1 ∂u ¯02 ∂s 2 ∂u ¯1 ∂u ¯0 ∂s 1 ∂u + + +2 2 = Re ∂y ∂x ∂x Re ∂x ∂y ∂y      3 η2 ∂ ∂u ¯01 ∂u ¯02 ∂ 2 h ∂u ¯01 ∂u ¯02 ∂ h + 2 + + + 6Re ∂t ∂y ∂x ∂t∂x ∂y ∂x ∂t2 ∂x  0    3  ¯1 ∂u ¯02 ∂ 2 h ∂u ¯01 ∂u ¯02 ∂ h ∂ ∂u +2 +2 +2 +4 ∂t ∂x ∂y ∂t∂y ∂x ∂y ∂t2 ∂y

939

(64)

(65)

(66) (67) (68) (69)

(70)

(71) (72) (73)

(74)

1 + O(η 4 ) + f¯W at z = s 2

(75)

0 0 0 1 f¯W = f¯W = f¯W = f¯W = 0 at z = s 1 2 3 3

(76)

940

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

4.2. First order approximation We can now consider a first order approximation: u ˜εi = u ¯0i + ε¯ u1i

(i = 1, 2),

u ˜ε3 = u ¯03 + ε¯ u13 + ε2 u ¯23 ,

p˜ε = p¯0 + ε¯ p1

(77)

0 The terms u ¯k3 (k = 0, 1, 2), p¯k (k = 0, 1) and σ ¯i3 (i = 1, 2) are known (see (57)–(61)), so the approximations of the vertical velocity and the pressure are:

    0 ∂u ¯02 ∂b ∂b ∂u ¯1 + +u ¯01 +u ¯02 u ˜ε3 = ε¯ u13 + ε2 u ¯23 = ε (b − z) ∂x ∂y ∂x ∂y   1   1 ∂u ¯1 ∂u ¯2 ∂b ∂b + +u ¯11 +u ¯12 + ε2 (b − z) ∂x ∂y ∂x ∂y     ε ε ∂u ˜2 ∂b ∂b ∂u ˜1 + +u ˜ε1 +u ˜ε2 = ε (b − z) ∂x ∂y ∂x ∂y  0    1  ¯1 ∂u ¯02 ¯1 ∂u ¯12 2 ∂u 2 ∂u p˜ε = p¯s − + + ε (s − z)(G − 2φ(cos ϕ)¯ u01 ) − + Re ∂x ∂y Re ∂x ∂y  ε  ε ∂u ˜2 2 ∂u ˜1 + = p¯s + ε(s − z)(G − 2φ(cos ϕ)¯ u01 ) − Re ∂x ∂y  ε  ˜1 ∂u ˜ε2 2 ∂u + + O(ε2 ) = p¯s + ε(s − z)(G − 2φ(cos ϕ)˜ uε1 ) − Re ∂x ∂y

(78)

(79) (80)

To obtain closed equations for u ¯01 and u ¯02 we have to get rid of the terms u ¯21 and u ¯22 in Eqs. (63) and (65). We can accomplish this with the following procedure. First integrate those equations in the vertical variable z from b to s:  2 0  ∂ p¯s 1 ∂ u ∂2u ∂2u ∂u ¯01 ¯01 ¯01 ¯1 ¯01 ¯02 0 ∂u 0 ∂u +u ¯1 +u ¯2 + − 3 − 2φ (sin ϕ) u ¯02 + +2 ∂t ∂x ∂y ∂x Re ∂x2 ∂y 2 ∂x∂y    2  0 2 0  η2 ∂ u ∂u ¯02 ∂ 2 u 1 ∂u ¯21  ¯1 ¯01 ∂u ¯1  ¯1 ∂ u + + = O(η 4 ) − − (81) hRe ∂z z=s ∂z z=b 3 ∂t ∂t∂x ∂t ∂t∂y   ∂ p¯s 1 ∂2u ∂2u ∂2u ∂u ¯02 ∂u ¯0 ∂u ¯0 ¯02 ¯02 ¯01 +u ¯01 2 + u ¯02 2 + − + 2φ (sin ϕ) u ¯01 + 3 + 2 ∂t ∂x ∂y ∂y Re ∂x2 ∂y 2 ∂y∂x    2 2  0 2 0  η2 ∂ u ∂u ¯02 ∂ 2 u ∂ u ¯1 ∂ u 1 ¯2  ∂2u ¯22  ¯2 ¯02 + + = O(η 4 ) − − (82) hRe ∂z 2 z=s ∂z 2 z=b 3 ∂t ∂t∂x ∂t ∂t∂y ∂u ¯2i at the top surface and at the bottom. Using (70) and (58), we write: ∂z   0   ∂u ¯02 ∂ ∂u ¯2i ∂b ∂b ∂u ¯1 1 (b − z) = Re σ ¯i3 + +u ¯02 +u ¯01 i = 1, 2 (83) − ∂z ∂xi ∂x ∂y ∂x ∂y

Then, we can find the value of

Next, substituting z by s and b, we have:    0    ! ∂u ¯1 ∂u ¯02 ∂u ¯2i  ∂ ∂u ¯2i  1  1  + − = Re σ ¯i3 z=s − σ ¯i3 z=b + h ∂z z=s ∂z z=b ∂xi ∂x ∂y

i = 1, 2

(84)

1 To obtain the value of σ ¯i3 at the top surface and at the bottom, boundary conditions (72), (73), (74) and (75) are used

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

    0   ∂u ¯02 ∂h ∂u ¯02 ∂h ∂u ¯21  ∂u ¯01 ∂u ¯21  ∂u ¯1 + + + − =2 2 ∂z z=s ∂z z=b ∂x ∂y ∂x ∂y ∂x ∂y    2   3 2 0 0 ∂ ∂u ¯1 ∂u ¯2 ∂ h ∂u ¯01 ∂u ¯02 η ∂ h 4 2 + +2 2 + + 6 ∂t ∂x ∂y ∂t∂x ∂x ∂y ∂t2 ∂x  0     ¯1 ∂u ¯02 ∂ 2 h ∂u ¯01 ∂u ¯02 ∂3h ∂ ∂u + + + +2 ∂t ∂y ∂x ∂t∂y ∂y ∂x ∂t2 ∂y  0  ! ¯1 ∂u ¯02 1 ¯1 + h ∂ ∂ u + + O(η 4 ) + f + Re f¯W R 1 1 ∂x ∂x ∂y     0   0 ∂u ¯02 ∂h ∂u ¯02 ∂h ∂u ¯22  ∂u ¯22  ∂u ¯1 ∂u ¯1 +2 + + − =2 ∂z z=s ∂z z=b ∂x ∂y ∂y ∂y ∂x ∂x     0  3 ∂ ∂u ∂u ¯0 ∂ 2 h ∂u ¯0 η2 ¯01 ∂u ¯1 ∂ h 4 +2 2 +2 +2 2 + 6 ∂t ∂x ∂y ∂t∂y ∂x ∂y ∂t2 ∂y  0   0  3  ∂ ∂u ¯1 ∂u ¯02 ∂ 2 h ∂u ¯1 ∂u ¯0 ∂ h +2 + + + 2 ∂t ∂y ∂x ∂t∂x ∂y ∂x ∂t2 ∂x  0  ! ∂u ¯02 ¯1 1 ¯1 + h ∂ ∂ u + + O(η 4 ) + f + Re f¯W R2 2 ∂y ∂x ∂y

941

(85)

(86)

When we substitute expressions (85)–(86) into (81)–(82), we obtain  0 2 0  η2 ∂ u ∂u ¯02 ∂ 2 u ∂u ¯01 ∂u ¯0 ∂u ¯0 ¯1 ¯01 ¯1 ∂ u +u ¯01 1 + u ¯02 1 + + ∂t ∂x ∂y 3 ∂t ∂t∂x ∂t ∂t∂y  2 0  1 ∂ u ∂2u ∂2u ¯1 ¯01 ¯02 ∂ p¯s + 4 + 2φ (sin ϕ) u ¯02 + +3 =− 2 2 ∂x Re ∂x ∂y ∂x∂y "    0  ∂u ¯0 ∂u ¯02 ∂h ∂u ¯1 ∂u ¯02 ∂h 1 2 2 1+ + + + hRe ∂x ∂y ∂x ∂y ∂x ∂y    2   3 2 0 0 ∂ ∂u ¯1 ∂u ¯2 ∂ h ∂u ¯01 ∂u ¯02 η ∂ h 4 2 + +2 2 + + 6 ∂t ∂x ∂y ∂t∂x ∂x ∂y ∂t2 ∂x  0    # ∂u ¯02 ∂ 2 h ∂u ¯02 ∂ ∂u ¯1 ∂u ¯01 ∂3h + + + +2 ∂t ∂y ∂x ∂t∂y ∂y ∂x ∂t2 ∂y ! 1 ¯1 fW1 + f¯R1 1 + O(η 4 ) + h  0 2 0  0 η2 ∂ u ∂u ¯02 ∂ 2 u ∂u ¯0 ∂u ¯0 ¯2 ¯02 ∂u ¯2 ¯1 ∂ u +u ¯01 2 + u ¯02 2 + + ∂t ∂x ∂y 3 ∂t ∂t∂x ∂t ∂t∂y  2 0  1 ∂ u ∂2u ∂2u ¯2 ¯0 ¯01 ∂ p¯s + − 2φ (sin ϕ) u ¯01 + 4 22 + 3 =− 2 ∂y Re ∂x ∂y ∂y∂x " 0   0  ∂u ¯1 ∂u ¯02 ∂h ∂u ¯1 ∂u ¯0 ∂h 1 + +2 +2 2 + hRe ∂y ∂x ∂x ∂x ∂y ∂y   0  2  0  3 2 0 ∂ ∂u ∂u ¯2 ∂ h ∂u ¯02 η ¯1 ∂u ¯1 ∂ h 2 + + + + 6 ∂t ∂y ∂x ∂t∂x ∂y ∂x ∂t2 ∂x  0    3 # ∂u ¯02 ∂ 2 h ∂u ¯02 ∂ ∂u ¯1 ∂u ¯01 ∂ h +2 +2 +2 +4 ∂t ∂x ∂y ∂t∂y ∂x ∂y ∂t2 ∂y ! 1 ¯1 fW2 + f¯R1 2 + O(η 4 ) + h

(87)

(88)

942

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

We repeat the process from Eqs. (64) and (66) to obtain closed equations for u ¯1i (i = 1, 2). We yield: ∂u ¯11 ∂u ¯1 ∂u ¯0 ∂u ¯1 ∂u ¯0 +u ¯01 1 + u ¯11 1 + u ¯02 1 + u ¯12 1 ∂t ∂x ∂x ∂y ∂y  0 2 1  2 1 2 0 η ∂u ¯1 ∂ u ∂u ¯1 ∂ u ∂u ¯02 ∂ 2 u ∂u ¯12 ∂ 2 u ¯1 ¯1 ¯11 ¯01 + + + + 3 ∂t ∂t∂x ∂t ∂t∂x ∂t ∂t∂y ∂t ∂t∂y   ∂s h ∂u ¯02 ∂h 0 ∂u ¯01 1 0 ∂b = − G + 2φ (sin ϕ) u ¯2 + 2φ(cos ϕ) u ¯ +h + −u ¯2 ∂x ∂x 1 ∂x 2 ∂y ∂y  2 1    ∂ u 2 ∂u ¯1 ∂u ¯12 ∂h ∂2u ∂2u ¯1 ¯11 ¯12 1 4 + 2 1+ + +3 + 2 2 Re ∂x ∂y ∂x∂y hRe ∂x ∂y ∂x  1   1 2 1 ∂u ¯2 ∂h η ∂h ∂2h ∂u ¯1 + + + 2 (2φ(cos ϕ)¯ u01 − G) hRe ∂y ∂x ∂y 6h ∂t ∂t∂x     4 ∂ ∂u ¯1 ∂u ¯12 ∂ 2 h 2 ∂u ¯1 ∂u ¯12 ∂3h + 2 1+ + 2 1+ Re ∂t ∂x ∂y ∂t∂x Re ∂x ∂y ∂t2 ∂x  1    3  2 ∂ ∂u ¯1 ∂u ¯12 ∂ 2 h 1 ∂u ¯11 ∂u ¯12 ∂ h + + + + Re ∂t ∂y ∂x ∂t∂y Re ∂y ∂x ∂t2 ∂y ! 1 ¯2 fW1 + f¯R2 1 + O(η 4 ) + h ∂u ¯1 ∂u ¯0 ∂u ¯1 ∂u ¯0 ∂u ¯12 +u ¯01 2 + u ¯11 2 + u ¯02 2 + u ¯12 2 ∂t ∂x ∂x ∂y ∂y  0 2 1  2 1 2 0 η ∂u ¯1 ∂ u ∂u ¯02 ∂ 2 u ∂u ¯12 ∂ 2 u ¯2 ¯2 ¯12 ¯02 ∂u ¯1 ∂ u + + + + 3 ∂t ∂t∂x ∂t ∂t∂x ∂t ∂t∂y ∂t ∂t∂y   ¯01 ∂s ∂s 0 h ∂ u ¯11 + 2φ(cos ϕ) u ¯1 + = − G − 2φ (sin ϕ) u ∂y ∂y 2 ∂y  2 1   1  2 1 2 1 2 ∂u ¯1 ∂u ¯12 ∂h ∂ u ∂ u ¯2 ¯2 ¯1 1 ∂ u + +2 +4 2 +3 + Re ∂x2 ∂y ∂x∂y hRe ∂x ∂y ∂y  1   1 ∂u ¯1 ∂u ¯12 ∂h η2 ∂h ∂2h + + + 2 (2φ(cos ϕ)¯ u01 − G) hRe ∂y ∂x ∂x 6h ∂t ∂t∂y  1  2  1  3 1 1 ∂u ¯ 2 ∂u ∂u ¯ 4 ∂ ∂u ¯1 ∂ h ¯1 ∂ h +2 2 + +2 2 + Re ∂t ∂x ∂y ∂t∂y Re ∂x ∂y ∂t2 ∂y  1    3  2 ∂ ∂u ∂u ¯12 ∂ 2 h 1 ∂u ∂u ¯12 ¯1 ¯11 ∂ h + + + + Re ∂t ∂y ∂x ∂t∂x Re ∂y ∂x ∂t2 ∂x ! 1 ¯2 f + f¯R2 2 + O(η 4 ) + h W2 Eqs. (87)–(90) are used to obtain the following equations for u ˜εi (see (77)):  ∂u ˜ε1 ∂ 2 u ∂u ˜ε2 ∂ 2 u ˜ε1 ˜ε1 + ∂t ∂t∂x ∂t ∂t∂y  2 ε  ∂s ∂ u 1 ∂2u ∂2u ˜1 ˜ε1 ˜ε2 ∂ p¯s ε − ε G + 2φ (sin ϕ) u ˜2 + 4 + +3 =− ∂x ∂x Re ∂x2 ∂y 2 ∂x∂y "    ε  ε ε ε ∂u ˜ ∂u ˜2 ∂h ∂u ˜1 ∂u ˜2 ∂h 1 2 2 1+ + + + hRe ∂x ∂y ∂x ∂y ∂x ∂y      3 ∂ ∂u ˜ε1 ∂u ˜ε2 ∂ 2 h ∂u ˜ε1 ∂u ˜ε2 η2 ∂ h 4 2 + +2 2 + + 6 ∂t ∂x ∂y ∂t∂x ∂x ∂y ∂t2 ∂x

η2 ∂u ˜ε1 ∂u ˜ε ∂u ˜ε +u ˜ε1 1 + u ˜ε2 1 + ∂t ∂x ∂y 3



(89)

(90)

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

# ∂3h ∂t2 ∂y   ∂h ε h ∂u ˜ε2 η 2 ∂h ∂ 2 h ∂u ˜ε1 ε ε ∂b + 2φ(cos ϕ)ε u ˜ +h + −u ˜2 + (2φ(cos ϕ)˜ u1 − G)ε 3h ∂t ∂t∂x ∂x 1 ∂x 2 ∂y ∂y ! 1 ˜ε fW1 + f˜Rε 1 + O(η 4 ) + O(ε2 ) + εh  ε 2 ε  ε η2 ∂ u ∂u ˜ε2 ∂ 2 u ∂u ˜ε ∂u ˜ε ˜2 ˜ε2 ∂u ˜2 ˜1 ∂ u +u ˜ε1 2 + u ˜ε2 2 + + ∂t ∂x ∂y 3 ∂t ∂t∂x ∂t ∂t∂y  2 ε  ∂s 1 ∂ u ∂2u ∂2u ˜2 ˜ε2 ˜ε1 ∂ p¯s ε − ε G − 2φ (sin ϕ) u ˜1 + +4 2 +3 =− ∂y ∂y Re ∂x2 ∂y ∂y∂x " ε   ε  ε ε ∂u ˜1 ∂u ˜2 ∂h ∂u ˜1 ∂u ˜ ∂h 1 + +2 +2 2 + hRe ∂y ∂x ∂x ∂x ∂y ∂y   ε   ε  3 η2 ∂ ∂u ˜1 ∂u ˜ε2 ∂ 2 h ∂u ˜1 ∂u ˜ε2 ∂ h + 2 + + + 6 ∂t ∂y ∂x ∂t∂x ∂y ∂x ∂t2 ∂x  ε   ε  3 # ∂u ˜ε ∂ 2 h ∂u ˜ε ∂ ∂u ˜1 ∂u ˜1 ∂ h +2 2 +2 +2 2 +4 ∂t ∂x ∂y ∂t∂y ∂x ∂y ∂t2 ∂y   ∂s ε ∂ u η 2 ∂h ∂ 2 h ˜ε1 h ε + 2φ(cos ϕ)ε u ˜ + + (2φ(cos ϕ)˜ u1 − G)ε 3h ∂t ∂t∂y ∂y 1 ∂y 2 ! 1 ˜ε f + f˜Rε 2 + O(ε2 ) + O(η 4 ) + εh W2 ∂ +2 ∂t



∂u ˜ε2 ∂u ˜ε1 + ∂y ∂x



∂2h + ∂t∂y



∂u ˜ε2 ∂u ˜ε1 + ∂y ∂x

943



(91)

(92)

The system of equations for u ˜i (i = 1, 2) (91)–(92) is coupled with the following equation for the water depth deduced from (67)–(68):   2 ε  2 ε  ∂ uε1 ) ∂(h˜ uε2 ) η 2 ∂ ∂ u ∂ u ∂h ∂(h˜ ˜ ˜ + + − h 21 + h 22 = O(η 4 ) ∂t ∂x ∂y 6 ∂x ∂t ∂y ∂t

(93)

Therefore, if Eqs. (91)–(93) are solved, we obtain the horizontal velocity and the free surface elevation. 5. Proposed shallow water model In this section, we go back the original variables: H ε (tε , xε , y ε ) = εLC h(t, x, y), u ˜εi (t, x, y) =

TC ˜ ε ε ε ε U (t , x , y ), LC i

p˜ε (t, x, y, z) =

B ε (xε , y ε ) = εLC b(x, y),

(94)

i = 1, 2, 3,

(95)

TC2 ˜ ε ε ε ε ε P (t , x , y , z ), ρ0 L2C

2 ˜f ε = TC F ˜ ε , W ρ0 L2C W

2 ˜f ε = TC F ˜ ε R ρ0 L2C R

p¯s (t, x, y) =

TC2 ¯ ε ε ε ε P (t , x , y ) ρ0 L2C s

(96) (97)

and we present the model that we have derived. We begin summarizing the results that we have achieved in this theorem: Theorem 1. Let us suppose that there exists asymptotic expansion (47). Then approximated solution (94)–(96) verifies

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

944

   ˜ ε ε ∂H ε η 2 TC2 ∂ 2 U εLC  ε ε ˜ +∇ · H U − = O(η 4 ) ε ε 2 ∂t 6 ∂(t ) TC ε

˜ 2 2 ∂U ˜ ε · U ˜ ε + η TC ∇ + ∇ U ∂tε 3



ε

˜ ∂U ∂tε

(98)

ε

˜   ∂U 1 ¯ε ˜ ε + 3ν∇ ∇ · U ˜ ε · = − ∇ P + νΔ U s ∂tε ρ0

  ˜ε   2 ε # ∂R η 2 TC2 ∂H ε ∂ H ε ε ε ˜ ˜ R · ∇H + +R ·∇ 2 ε ·∇ ε 6 ∂t ∂t ∂(tε )2  ˜ε     ! η 2 TC2 ∂H ε ∂H ε U2 ε ε ˜ + 2Φ (sin ϕ) ˜ ε + 2Φ(cos ϕ)U1 − g ∇S + 3H ε ∂tε ∇ ∂tε −U 1   ε ˜ ε  ˜ · ∇B ε + 1 H ε ∇ · U 1 −U + 2Φ(cos ϕ) H ε ∇(U˜1ε ) + 2 2 0  1  ˜ ε ˜ ε + LC O(η 4 ) + LC O(ε2 ) + FW + F R ε ρ0 H TC2 TC2 ν + ε H

"

(99)

2 ! ˜ ε + O(ε2 ) ρ0 LC ˜ ε − 2μ∇ε · U P˜ ε = P¯s + ρ0 (S ε − z ε ) g − 2Φ(cos ϕ)U 1 2 TC

(100)

˜ ε + U ˜ ε · ∇B ε ˜ ε = (B ε − z ε ) ∇ε · U U 3

(101)

˜ ε = (U ˜ ε, U ˜ ε ) is the time-averaged horizontal velocity approximation, P¯s is the time-averaged atmowhere U 1 2 ˜ ε , F ˜ ε are the wind and friction forces and spheric pressure at the surface, F W

R

⎛

˜ε ˜ε ∂U ∂U 1 2 + 2 ε ⎜ ∂xε ∂y ε ˜ =⎜ R ⎝ ∂U ˜ε ˜ε ∂U 1 2 + ε ∂y ∂xε 4

˜ε ˜ε ∂U ∂U 1 2 + ∂y ε ∂xε ˜ε ˜ε ∂U ∂U 2 1ε + 4 ε2 ∂x ∂y

⎞ ⎟ ⎟ ⎠

We, finally, propose the following model that we have derived neglecting the O(ε2 ) and O(η 4 ) terms in the above equations. For notational convenience, we henceforth drop the ˜.    ¯ ε ε ∂H ε η 2 TC2 ∂ 2 U  ε ε ¯ +∇ · H U − =0 ∂tε 6 ∂(tε )2 ε

¯ 2 2 ∂U ¯ ε · U ¯ ε + η TC ∇ + ∇ U ∂tε 3



ε

¯ ∂U ∂tε

ε

¯   ∂U 1 ¯ε ¯ ε + 3ν∇ ∇ · U ¯ ε · = − ∇ P + νΔ U s ∂tε ρ0

  2 ε #  ¯ε  η 2 TC2 ∂H ε ∂ H ∂R ε ε ε ¯ ¯ R · ∇H + +R ·∇ 2 ε ·∇ ε 6 ∂t ∂t ∂(tε )2  ¯ε     ! η 2 TC2 ∂H ε ∂H ε U2 ε ε ¯ + 2Φ (sin ϕ) ¯ ε + 2Φ(cos ϕ)U1 − g ∇S + 3H ε ∂tε ∇ ∂tε −U 1   ε  ¯ ε  ¯ · ∇B ε + 1 H ε ∇ · U 1 1  ¯ ε −U ¯ ε + 2Φ(cos ϕ) H ε ∇(U¯1ε ) + + F + F 2 W R 2 ρ0 H ε 0 ! ¯ ε ¯ ε − 2μ∇ε · U P¯ ε = P¯s + ρ0 (S ε − z ε ) g − 2Φ(cos ϕ)U 1 ν + ε H

(102)

"

¯ ε + U ¯ ε · ∇B ε ¯ ε = (B ε − z ε ) ∇ε · U U 3

(103) (104) (105)

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

945

Remark 1. In practice, we neglect the term in η 2 in Eq. (102), because it does not provide significant improvements in the accuracy of the numerical solutions that we present in Section 6, but it causes some problems of stability. Therefore, in what follows, we replace (102) by   ∂H ε ¯ ε = 0 + ∇ε · H ε U ε ∂t

(106)

6. Numerical experiments In this section we shall compare model (103)–(106) (we shall refer to it as NM) with other shallow water model, not averaged in time, obtained in [8] using the same asymptotic techniques. We shall refer to it as SW, and it can be written as in [9]:   ∂H ε ε ε ε + ∇ · H =0 U ∂tε ε ∂U ε·U  ε = − 1 ∇P ε + ν{ΔU  ε + 1 [∇U  ε T + ∇U  ε ]∇H ε + ∇U s ε ε ∂t ρ0 H "  ε    1 ε 2  ε )]} − g∇S ε + 1  ε + 2Φ (sin ϕ) U2 ε +F + ∇[(H ) (∇ · U F W R −U1ε (H ε )2 ρ0 H ε #   ε  ε · ∇B ε + 1 H ε ∇ · U 1 −U (107) + (cos ϕ) U1 ∇S ε + H ε ∇(U1ε ) + 2 2 0  ε is not averaged in time. where the horizontal velocity U In order to compare models (103)–(106) and (107), we shall consider some analytical solutions of Navier–Stokes equations (1)–(8) whose velocity rapidly oscillates in time. Then we solve numerically Eqs. (103)–(106) and (107) for the data provided by the analytical solutions of Navier–Stokes equations, and finally we compute the errors committed by each of the models. To perform the numerical simulations, we have opted to use MacCormack scheme (see [6]) due to its good stability properties, its easy implementation, and to the fact that has been applied successfully to the resolution of similar problems. Let us introduce now the first family of exact solutions to Navier–Stokes equations that we shall use to compare models SW and NM. Horizontal and vertical velocity oscillate in time, but the horizontal velocity depends on variable xε while the vertical velocity depends on variable z ε :    2πn1 ε 2πn1 ε t + (B1 + B2 xε ) cos t Tp Tp     2πn2 ε 2πn2 ε ε ε t + (D1 + D2 x ) cos t + (C1 + C2 x ) sin Tp Tp 

U1ε = (A1 + A2 xε ) sin

U2ε = 0

         2πn1 ε 2πn1 ε 2πn2 ε 2πn2 ε t + B2 cos t + C2 sin t + D2 cos t U3ε = −z ε A2 sin Tp Tp Tp Tp

Bε = 0 H

ε

= Ee

Tp 2π

$

A2 n1

       % 2πn B 2πn C 2πn D 2πn cos Tp 1 tε − n 2 sin Tp 1 tε + n 2 cos Tp 2 tε − n 2 sin Tp 2 tε

Psε = P ε (z ε = H ε ) − 2μ

1

∂U3ε ε + FW 3 ∂z ε

2

2

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

946

Table 1 Error bounds for example (108) with data (109) and ηTC = 2.5. Δt

Error bound for H SW

Error bound for H NM

Error bound for U1 SW

Error bound ¯1 NM for U

0.01 0.005 0.001 0.0001 0.000025

0 0 0 0 0

0 0 0 0 0

3.1e3 2.6e0 4.0e−1 3.8e−2 9.5e−3

9.7e−3 4.8e−3 9.7e−4 9.7e−5 2.4e−5

ε ε FW = FW =0 1 2

FRε = 0

(108)

where Ai , Bi , Ci , Di , ni (i = 1, 2) and Tp are any real value. We are able to calculate an analytical expression for P ε from Eq. (1), but it is too long and we have decided not to include it here. We consider that D is a rectangular basin of length 10 meters and width 2 meters with a 100 × 20 points grid (that is, the discretization step used is Δxε = Δy ε = 0.1). We choose the values of the parameters so the maximum depth is always smaller than 1 meter, and thus the aspect ratio is always smaller than 10−1 . We solve in temporal interval [0, 10] with different time steps. We introduce these two sets of values for the constants: A1 = B1 = 2, C1 = D1 = 0.5, A2 = B2 = C2 = D2 = 0, E = 1

(109)

A1 = B1 = 1, C1 = D1 = 0.5, A2 = B2 = 0.5, C2 = D2 = 0, E = 0.75

(110)

With election (109) U1ε does not depend on xε , but with election (110) it does. We present in Tables 1 and 3 the infinity-norm errors obtained when we approximate solution (108) using models (103)–(106) and (107) with these choices of constants. For both examples we take: Φ = 0, g = 9.8, ρ = 998.2, ν = 1.02 × 10−6 , Tp = 1, n1 = 0.1, n2 = 100. We remark that in model (103)–(106) we must choose η and TC , where TC is the characteristic time and 2η is the length of the averaging interval for the non-dimensional time variable; but the relevant choice here is the value of the product ηTC because 2ηTC is the length of the averaging interval for the dimensional time variable. In Table 1 we show the error bounds when comparing the computed solution of SW with the Navier– Stokes solution (108), and when comparing the computed solution of NM with the averaged version of (108). We can see that, to obtain the same order of accuracy than the new model (103)–(106) with Δt = 10−2 , model (107) has to be solved taking Δt = 2.5 × 10−5 , that is, the time step must be 400 times smaller. In this example, the depth is computed exactly by both models because it just depends on time. The times (in seconds) required to solve models SW and NM for the example of Table 1 are presented in Table 2. We have used for these computations a personal computer with an Intel(R) Core(TM) i5 2.80 GHz processor (6 GB RAM). We can see in this Table 2 that execution time for the new model increased compared with model SW for the same time step, but we can observe that it is shorter if we want to obtain the same accuracy. For example in this case, using the new model, in 8.4 seconds we have a more accurate approximation of the solution than the one obtained with model SW in 268 seconds. We see that the time needed to achieve the same accuracy with (107) is 1073 seconds. In Table 3 we observe that, when velocity depends on the spatial variable xε , the results achieved with the new model are quite better too. When we introduce data (110), we need to introduce a small value for η (ηTC = 5 × 10−3 ) to guarantee the convergence of the scheme for Δt = 10−4 . For larger values of Δt, it is possible to choose larger values of η.

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

947

Table 2 Execution times (in seconds) of the example of Table 1.

SW NM

Δt = 10−2

Δt = 5 × 10−3

Δt = 10−3

Δt = 10−4

Δt = 2.5 × 10−5

2.9 8.4

5.4 16.4

27.0 81.5

268.5 812.6

1073.0 3284.9

Table 3 Error bounds for example (108) with data (110) and ηTC = 5 × 10−3 . Δt

Error bound for H SW

Error bound for H NM

Error bound for U1 SW

¯1 NM Error bound for U

0.005 0.001 0.0001

9.8e−3 1.9e−3 1.9e−4

9.8e−3 1.9e−3 1.9e−4

3.1e0 4.7e−1 4.4e−2

1.4e−2 2.7e−3 2.7e−4

Fig. 2. First order accuracy of the numerical scheme.

In Fig. 2 we show the error bounds of Table 3. It can be well appreciated that the numerical scheme is first order accurate as it corresponds to MacCormack scheme. Let us consider now the following solution to Navier–Stokes equations where water depth depends on xε : U1ε

         2πn1 ε 2πn1 ε 2πn2 ε 2πn2 ε = Ax sin t + cos t + sin t + cos t Tp Tp Tp Tp ε

U2ε = 0

         2πn1 ε 2πn1 ε 2πn2 ε 2πn2 ε U3ε = −Az ε sin t + cos t + sin t + cos t Tp Tp Tp Tp Bε = 0

⎡

H ε = B(xε )m e

ATp (m+1) ⎣ 2π

Psε = P ε (z ε = H ε ) − 2μ

      ⎤  2πn1 ε 2πn2 ε 2πn1 ε 2πn2 ε cos cos t t t t −sin −sin Tp Tp Tp Tc ⎦ + n1 n2

∂U3ε ε + FW 3 ∂z ε

ε ε FW = FW =0 1 2

FRε = 0

(111)

with m ∈ N, A, B, ni (i = 1, 2) and Tp any real value. For the following values of the constants A = 0.03, B = 0.1, n1 = 0.1, n2 = 100, Tp = 1, m = 1

(112)

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

948

Table 4 Error bounds for example (111) with data (112) and ηTC = 1.3034 × 10−2 . Δt

Error bound for H SW

Error bound for H NM

Error bound for U1 SW

¯1 NM Error bound for U

0.001 0.0001

– 7.2e−2

1.2e−2 2.9e−3

– 3.6e−1

2.1e−1 5.1e−2

¯1 NM. Fig. 3. Comparison between U1 SW and U

time step must be very small (Δt = 10−4 ) to achieve the convergence of model SW while the new model gives reasonable results with Δt = 10−3 as we show in Table 4. We can observe from Tables 3 and 4 that, when the solution depends on xε , we need to take small values of ηTC . This is due to the fact that approximation (31) (and, as consequence, all formula (32)–(36)) is only valid for small values of η (specially when the approximated function is rapidly oscillating), but then, averaged solution is almost equal to exact solution, and we shall also need a very small time step for model (103)–(106). Nevertheless, we see in Table 4 that model (103)–(106) achieves better results than the model (107) even in the worst case. Finally, we present a more realistic numerical experiment which does not give a precise solution of the Navier–Stokes equations. Zero–Dirichlet boundary conditions are imposed for horizontal velocity. At the initial time, horizontal velocity is set to be zero and the water depth is:  H0ε (x, y) =

1 + 0.1 sin(2πx)

if x ∈ [0, 3) ∪ (7, 10]

1.9 + 0.1 sin(2πx) if [3, 7]

(113)

In this case, as we have already commented for the previous example, time step must be quite small (Δt = 5 × 10−4 ) to achieve the convergence of model SW while the new model gives reasonable results even with Δt = 10−2 . We have plotted together, in Fig. 3, the approximations of the U1 component of the horizontal velocity that both models provide. We observe that the new model reduces the high frequency oscillations that appear using the shallow water model. If we compare the execution times, we have that the time required to solve model SW with this data is 45 seconds while model NM runs in just 7 seconds. 7. Conclusions We have used asymptotic analysis to obtain from the time-averaged non-dimensional Navier–Stokes equations a new shallow water model. The new model is able to filter (in some representative cases) the high frequency oscillations and this allows us to choose a much larger time step. We have made some numerical comparisons between the new model and the shallow water model proposed in [8]. Numerical experiments confirm that this new model is able to obtain a given accuracy using larger time steps than the time step needed by the other shallow water model (see Tables 1, 3–4 and Fig. 3). In some cases, the time step can be even four hundred times larger. This enhancement leads to much shorter

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

949

execution times compared with the execution times required by the classical shallow water model to obtain the same precision (see Table 2). Although the numerical results achieved improve those of the model without time filtering, we have observed not so good results in some cases. We think that it can be caused by the absence of spacial filtering, and it would be convenient to use a combination of spacial and time filters, as suggested in [2], because the use of spacial “projective” filters or modeling of subgride-scale stress tensor is necessary to reduce the number of degrees of freedom of the problem. References [1] D. Carati, F.S. Winckelmans, H. Jeanmart, On the modelling of the subgrid-scale and filtered-scale stress tensors in large-eddy simulation, J. Fluid Mech. 441 (2001) 119–138. [2] D. Carati, A.A. Wray, Time filtering in large eddy simulations, in: Proceedings of the Summer Program 2000 Center for Turbulence Research, 2000, pp. 263–270. [3] M.J. Castro, J. Macías, Modelo matemático de las corrientes forzadas por el viento en el mar de Alborán, Grupo de Análisis Matemático, Universidad de Málaga, 1994. [4] J.-L. Guermond, J.T. Oden, S. Prudhomme, Mathematical perspectives on large eddy simulation models, J. Math. Fluid Mech. Turb. Flows 6 (2004) 194–248. [5] P.L. Lions, Mathematical Topics in Fluid Mechanics. Vol. 1: Incompressible Models, Oxford University Press, 1996. [6] R.W. MacCormack, Numerical solution of the interaction of a shock wave with a laminar boundary layer, in: M. Holt (Ed.), Proceedings 2nd Int. Conf. on Num. Methods in Fluid Dynamics, Springer-Verlag, 1971, pp. 151–163. [7] C. Pruett, Temporal large-eddy simulation: theory and implementation, Theor. Comput. Fluid Dyn. 22 (2007) 275–304. [8] J.M. Rodríguez, R. Taboada-Vázquez, From Navier–Stokes equations to shallow waters with viscosity by asymptotic analysis, Asymptot. Anal. 43 (4) (2005) 267–285. [9] J.M. Rodríguez, R. Taboada-Vázquez, Comportamiento numérico de un modelo de aguas someras con viscosidad, in: Proceedings of the Conference Métodos Numéricos en Ingeniería, 2005. [10] J.M. Rodríguez, R. Taboada-Vázquez, From Euler and Navier–Stokes equations to shallow waters with viscosity by asymptotic analysis, Adv. Eng. Softw. 38 (2007) 399–409. [11] J.M. Rodríguez, R. Taboada-Vázquez, A new shallow water model with linear dependence on depth, Math. Comput. Modelling 48 (2008) 634–655, http://dx.doi.org/10.1016/j.mcm.2007.11.002. [12] J.M. Rodríguez, R. Taboada-Vázquez, A new shallow water model with polynomial dependence on depth, Math. Methods Appl. Sci. 31 (2008) 529–549, http://dx.doi.org/10.1002/mma.924. [13] J.M. Rodríguez, R. Taboada-Vázquez, Bidimensional shallow water model with polynomial dependence on depth through vorticity, J. Math. Anal. Appl. 359 (2) (2009) 556–569. [14] J.M. Rodríguez, R. Taboada-Vázquez, Derivation of a new asymptotic viscous shallow water model with dependence on depth, Appl. Math. Comput. 219 (2012) 3292–3307. [15] P. Sagaut, Large Eddy Simulation for Incompressible Flows, Springer-Verlag, 2006. [16] R. Taboada-Vázquez, Modelos de aguas poco profundas obtenidos mediante la técnica de desarrollos asintóticos, PhD thesis, Universidade da Coruña, 2006. [17] Tan Weiyan, Shallow Water Hydrodynamics, Elsevier, 1992.

Further reading [18] E. Audusse, M.O. Bristeau, B. Perthame, Kinetic schemes for Saint-Venant equations with source terms on unstructured grids, Rapp. Rech. INRIA 3989 (2000) 1–44, https://hal.archives-ouvertes.fr/file/index/docid/72657/filename/RR-3989. ps. [19] P. Azérad, F. Guillén, Mathematical justification of the hydrostatic approximation in the primitive equations of geophysical fluid dynamics, SIAM J. Math. Anal. 33 (4) (2001) 847–859, http://dx.doi.org/10.1137/S0036141000375962. [20] E. Casas, Introducción a las ecuaciones en derivadas parciales, Universidad de Cantabria, ISBN 8487412750, 1992. [21] M.J. Castro, J.M. González-Vida, C. Parés, Numerical treatment of wet/dry fronts in shallow flows with modified roe schemes, Math. Models Methods Appl. Sci. 16 (6) (2006) 897–932, http://dx.doi.org/10.1142/S021820250600139X. [22] P.G. Ciarlet, Mathematical Elasticity. Volume II: Theory of Plates, North-Holland, ISBN 0444825703, 1997. [23] P.G. Ciarlet, Mathematical Elasticity. Volume III: Theory of Shells, North-Holland, ISBN 0444828915, 2000. [24] P. Destuynder, Une Théorie Asymptotique des Plaques Minces in Élasticité Linéaire, Masson, ISBN 2225807728, 1986. [25] J.-F. Gerbeau, B. Perthame, Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation, Discrete Contin. Dyn. Syst. Ser. B 1 (1) (2001) 89–102, http://dx.doi.org/10.3934/dcdsb.2001.1.89. [26] B. Di Martino, P. Orenga, M. Peybernes, Simulation of a spilled oil slick with a shallow water model with free boundary, Math. Models Methods Appl. Sci. 17 (3) (2007) 393–410, http://dx.doi.org/10.1142/S0218202507001966. [27] R. Teman, A. Miranville, Mathematical Modeling in Continuum Mechanics, Cambridge University Press, ISBN 0521643627, 2001.

950

J.M. Rodríguez, R. Taboada-Vázquez / J. Math. Anal. Appl. 428 (2015) 930–950

[28] L. Trabucho, J.M. Viaño, Mathematical modelling of rods, in: P.G. Ciarlet, J.-L. Lions (Eds.), Handbook of Numerical Analysis, vol. IV, North-Holland, 1996, pp. 487–974. [29] G.B. Whitham, Linear and Nonlinear Waves, John Wiley & Sons, ISBN 0471359424, 1974. [30] R.K. Zeytounian, Modélisation asymptotique en mécanique des fluides newtoniens, Springer-Verlag, ISBN 3540578382, 1994.