Hopf bifurcation in viscous incompressible flow down an ... - IUPUI Math
J. math. fluid. mech. 99 (9999), 1–22 1422-6928/99000-0, DOI 10.1007/s00009-003-0000 c 2006 Birkh¨
auser Verlag Basel/Switzerland
Journal of Mathematical Fluid Mechanics
Hopf bifurcation in viscous incompressible flow down an inclined plane: a numerical approach Roland Glowinski and Giovanna Guidoboni Abstract. In this article we investigate numerically a Hopf bifurcation phenomenon for a viscous incompressible flow down an inclined plane. This problem has been discussed by Nishida et al. who proved the existence of periodic solutions bifurcating from the steady flow. Using a computational methodology based on finite elements for the space discretization and on operator splitting for the time discretization, we have been able to reproduce the results predicted by Nishida et al. Mathematics Subject Classification (2000). 70K50, 76B45, 76M10. Keywords. free surface, surface tension, Hopf bifurcation, operator splitting, finite elements.
1. Introduction In [15], Nishida, Teramoto and Yoshihara proved the existence of a Hopf bifurcation for the flow of an incompressible viscous fluid down an inclined plane. The main goal of this article is to investigate numerically the same flow problem in order to precise at which value of the Reynolds number the Hopf bifurcation is occurring, and to visualize the shape of the past-bifurcation waves. The computational methodology discussed in this article relies on the following components: b (i) As in [9] we map the two dimensional flow region Ω(t) on a rectangular one Ω (actually the main result in [15] is proved using a similar transformation, which is quite classical in this context). b and we use a time discretization by (ii) We formulate the flow problem in Ω operator splitting to reduce the solution of the original flow problem to that of simpler sub-problems. Actually, it is easier to solve some of the sub-problems in
2
Roland Glowinski and Giovanna Guidoboni
JMFM
the physical flow region, and we take advantage of this simplification each time it is possible. (iii) We combine the time-splitting scheme with an isoparametric Bercoveir-Pironneau finite element approximation of the Navier-Stokes equations (see e.g. [3], [8], [9], [10], [18]) to fully discretize the free surface problem. (iv) The wave-like equation methodology discussed in, e.g., [8], [9], [10], and [17], is used to handle the pure advection problems resulting from the time-splitting and from the domain transformation. The resulting methodology is modular, relatively easy to implement, and it introduces very little numerical dissipation, an important property if one wishes to identify accurately the critical Reynolds number at which bifurcation occurs. We took advantage of this property in, e.g., [10] to study other Hopf bifurcation phenomena. The computed results are in very good agreement with the experimental data reported in [13]. The flow problem in this article is a basic one. It is not surprising therefore that it has attracted the attention of many investigators. Let us mention, among many others, [1], [2], [5], [6], [11] and [16]. The content of the remainder of this article is as follows. In Section 2 we describe the mathematical formulation of the problem. In Section 3 we discuss the time discretization by operator splitting. In Sections 4 and 5 the splitting scheme of Section 3 is combined with the Bercovier-Pironneau finite element approximation of the Navier-Stokes equations to derive the full discretization of the problem. The results of the numerical experiments are presented in Section 6.
2. Mathematical formulation An incompressible viscous fluid flows down an inclined plane under the action of gravity. Let us consider a Cartesian frame of reference with the x1 -axis (resp. x2 -axis) along (resp. perpendicular to) the plane which is inclined at an angle β, as shown in Figure 1. We assume that the free surface Γ(t) admits a cartesian representation at each instant of time so that it can be described by the timedependent function η(x1 , t). Therefore, the domain occupied by the fluid is Ω(t) = {x = (x1 , x2 ) ∈ R2 : x1 ∈ (0, L), x2 ∈ (0, η(x1 , t))},
(2.1)
the rigid lower boundary is Σ = {x = (x1 , x2 ) ∈ R2 : x1 ∈ (0, L), x2 = 0},
(2.2)
and the upper free surface is Γ(t) = {x = (x1 , x2 ) ∈ R2 : x1 ∈ (0, L), x2 = η(x1 , t)}.
(2.3)
The fluid flow is modeled by the Navier-Stokes equations ∇·u =0 ̺[∂t u + (u · ∇)u] = ∇ · σ + ̺g,
(2.4)
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
x
3
2
Γ(t)
η(x1,t)
0
β L
x
1
Figure 1. Sketch of the geometry of the problem. where u = (u1 , u2 ) is the velocity of the fluid, ̺ is the fluid density, σ is the stress tensor defined as σ = −pI + 2µD(u), where p is the pressure, µ is the viscosity and D(u) = (∇u + (∇u)t )/2 is the rate-of-strain tensor, where the superscript t means transpose. The external force is due to the gravity and it is given by ̺g = ̺g(sin β, − cos β), where g is the acceleration of gravity. Here ∂t denotes the partial derivative with respect to time, while ∇ = (∂1 , ∂2 ) is the gradient with respect to the spatial variables. Equations (2.4) hold in the time-dependent domain Ω(t) defined in (2.1), and for 0 < t < T . No-slip condition is assumed for the velocity at the rigid bottom, and therefore we have u=0
on Σ.
(2.5)
The free boundary Γ(t) is a material surface, therefore we prescribe the following kinematic condition ∂t η + u1 ∂1 η − u2 = 0 on Γ(t), (2.6) to ensure that the fluid particles do not cross the surface. An additional condition is needed to establish the balance of stress σn = sH(η(x1 , t))n on Γ(t), (2.7) p 2 where n = (−∂1 η, 1)/ 1 + |∂1 η| is the outward unit normal vector, s is the surface tension (which is assumed to be positive and constant), and H(η(x1 , t)) = ∂12 η/(1 + |∂1 η|2 )3/2 is the curvature. We assume periodicity in the x1 -direction of given period L. The incompressibility of the fluid and the periodic boundary conditions imply that the area of the fluid domain must be conserved. Indeed we have: Z Z L Z Z d L η dx1 . (2.8) ∂t η dx1 = u · n dΓ(t) = ∇ · ud x = 0= dt 0 0 Γ(t) Ω(t) This means that
Z
0
L
η(x1 , t) dx1 = constant = HL,
(2.9)
4
Roland Glowinski and Giovanna Guidoboni
JMFM
where H is the average height of the fluid layer. To complete the formulation of the problem, initial conditions should be prescribed for the velocity of the fluid and the shape of the domain: u(x, 0) = u0 (x)
with ∇ · u0 = 0, Z L η(x1 , 0) = η0 (x1 ) with η0 (x1 ) dx1 = HL.
(2.10)
0
In order to write the variational formulation of the problem (2.4)-(2.10) we introduce the following functional spaces: V0 (Ω(t))
=
L20 (Ω(t))
=
S
{v | v ∈ (H 1 (Ω(t)))2 , v = 0 on Σ, v|x1 =0 = v|x1 =L } , (2.11) Z qdx = 0} , (2.12) {q | q ∈ L2 (Ω(t)), Ω(t)
=
1
{η | η ∈ H (0, L), η|x1 =0 = η|x1 =L } .
(2.13)
Now, we multiply equations (2.4)1 and (2.4)2 by the test functions q ∈ L2 (Ω(t)) and v ∈ V0 (Ω(t)), respectively, and we use the notation ϕ(t) to denote the function x → ϕ(x, t). The resulting equations are integrated over Ω(t) to obtain the variational formulation of (2.4)-(2.10): For a.e. t > 0, find u(t) ∈ V0 (Ω(t)), p(t) ∈ L20 (Ω(t)), η(t) ∈ S such that Z Z Z ̺ ∂t u · v dx + ̺ (u · ∇)u · v dx + 2µ D(u) : D(v) dx Ω(t) Ω(t) Ω(t) Z Z − p∇ · v dx = ̺ g · v dx Ω(t) Ω(t) Z +s H(η(x1 , t))n · v dΓ(t), ∀v ∈ V0 (Ω(t)), Γ(t) Z q∇ · u dx = 0, ∀q ∈ L2 (Ω(t)), Ω(t) Z L Z L ∂t η φ dx1 + (u1 |Γ(t) ∂1 η − u2 |Γ(t) ) φ dx1 = 0, ∀φ ∈ S, 0
(2.14)
0
u(x, 0) = u0 (x), with ∇ · u0 = 0, Z L η(x1 , 0) = η0 (x1 ) with η0 (x1 ) dx1 = HL. 0
Remark 1. There is no contradiction between the choice of S given by (2.13) and the fact that H(η) includes ∂12 η. Indeed, by integration by parts we get Z Z L s H(η)v1 n1 dΓ(t) = −s (1 + |∂1 η|2 )−1/2 ∂1 v1 dx1 , (2.15) Γ(t)
0
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
5
At
H
0
L
0
Ω (t)
L
Ω
−1
At
Figure 2. At maps the the current domain Ω(t) into the referb ence domain Ω. s
Z
H(η)v2 n2 dΓ(t) = −s
Z
L
∂1 η(1 + |∂1 η|2 )−1/2 ∂1 v2 dx1 ,
(2.16)
0
Γ(t)
where the x1 -differentiation on η is transferred on the two components of the test function v. Problem (2.14) is defined on the domain Ω(t) which changes in time and it is not known a priori. Therefore we are going to reformulate problem (2.14) into a fixed reference domain. Let At be a family of mappings which at each time t ∈ (0, T ) maps the b = (0, L) × (0, H) defined by current domain Ω(t) into the reference domain Ω At :
Ω(t) ∈ R2
b ∈ R2 → Ω
x = (x1 , x2 ) → ξ = (ξ1 , ξ2 ) = At (x) =
ξ1 = x1 ξ2 =
H x2 , η(x1 , t)
(2.17)
as sketched in Figure 2. We observe that the free surface Γ(t) is mapped into ˆ = {ξ ∈ R2 | ξ1 ∈ (0, L), ξ2 = H}. Γ The Jacobian matrix J of this transformation 1 0 J = J(η) = ξ2 H − ∂1 η(ξ1 , t) η(ξ1 , t) η(ξ1 , t)
(2.18)
and its determinant J are given by ,
J = J(η) =
H . η(ξ1 , t)
(2.19)
It is clear that this transformation is well defined as long as η(ξ1 , t) > 0 or, in other words, as long as the free surface does not touch the bottom. Let f = f (x, t) be a function defined on Ω(t) and fˆ = fˆ(ξ, t) = f (A−1 t (ξ), t) the correspondent
6
Roland Glowinski and Giovanna Guidoboni
b Due to the change of variables, we have function defined on Ω. Z Z fˆ(ξ, t) Jdξ. f (x, t) dx = b Ω
Ω(t)
Moreover, it follows from the chain rule that b fˆ, ∂t f = ∂t fˆ − w · ∇
JMFM
(2.20)
(2.21)
b = J−t ∇ and the domain velocity w is given by where ∇ ξ2 w(ξ1 , t) = ∂t x = ∂t η(ξ1 , t)e2 H (2.22) ξ2 = u2 (ξ1 , η(ξ1 , t), t) − u1 (ξ1 , η(ξ1 , t), t)∂1 η(ξ1 , t) e2 . H The last step is a consequence of the kinematic condition (2.6) and we used the notation e2 = (0, 1). Using (2.20) and (2.21) , problem (2.14) becomes: b pˆ(t) ∈ L2 (Ω), b η(t) ∈ S, such that ˆ (t) ∈ V0 (Ω), For a.e. t > 0, find u 0
̺
Z b u·v b u·v ˆ J(η)dξ ˆ J(η)dξ + ̺ (ˆ ˆ·v ˆ J(η)dξ − ̺ (w · ∇)ˆ u · ∇)ˆ ∂t u b b b Ω Ω Ω Z Z Z ˆ ·v b u) : D(ˆ b v) J(η)dξ − ˆ·v ˆ J(η)dξ ˆ J(η)dξ = ̺ g pˆ∇ +2µ D(ˆ Z
Z
b Ω
Z
b Ω
b Ω
b H(η) vˆ2 (ξ1 , H, t) − vˆ1 (ξ1 , H, t)∂1 η dξ1 , ∀ˆ v ∈ V0 (Ω), +s 0 Z b ·u b ˆ J(η)dξ = 0, ∀ˆ qˆ∇ q ∈ L2 (Ω), b Ω Z L Z L (ˆ u1 (ξ1 , H, t)∂1 η − uˆ2 (ξ1 , H, t)) φ dξ1 = 0, ∀φ ∈ S, ∂t η φ dξ1 + 0
L
(2.23)
0
b ·u ˆ (ξ, 0) = u ˆ 0 (ξ), with ∇ ˆ 0 J(η0 ) = 0, u Z L η(ξ1 , 0) = η0 (ξ1 ) with η0 (ξ1 ) dξ1 = HL. 0
b u) = (∇ˆ b u + (∇ˆ b u)t )/2. Here we used the notation D(ˆ
3. Time discretization The numerical solution of problem (2.23) contains three main difficulties: (a) the incompressibility condition and the related unknown pressure, (b) an advection term associated to the velocity of the fluid, (c) an advection term associated to the velocity of the domain. The difficulty associated with the unknown location of the free boundary has been removed through the transformation defined by (2.17). A time discretization by operator splitting allows to decouple the difficulties of the problem and to treat each of them with appropriate algorithms. Moreover it
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
7
makes possible to update the location of the free boundary as a sub-step, avoiding costly sub-iterations between the location of the free surface and the solution of the Navier-Stokes equations. Among the many possible versions, we advocate a very simple operator splitting method which is first-order accurate [14]. The low order accuracy is compensated, actually, by easy implementation and less cost in computational time. ˆn = u ˆ (tn ), pˆn = pˆ(tn ), Let ∆t be the time discretization step and tn = n∆t, u n n n ˆ and η being known, the scheme consists of solving the η = η(t ). For n ≥ 0, u following problems. n
b and pˆn+1/5 ∈ L2 (Ω) b such that ˆ n+1/5 ∈ V0 (Ω) 1. Find u 0
Z Z n+1/5 ˆ ˆn u −u n b un+1/5 ) : D(ˆ b v)J(η n )dξ ˆ · vJ(η )dξ + µ D(ˆ ̺ ∆t b b Ω Ω Z Z n+1/5 b n ˆ·v ˆ J(η n )dξ ˆ g p ˆ ∇ · v J(η )dξ = ̺ − b b Ω Ω Z L n n b , +s dξ1 , ∀ˆ v ∈ V0 (Ω) H(η ) v ˆ (ξ , H, t) − v ˆ (ξ , H, t)∂ η 2 1 1 1 1 0 Z b ·u b . ˆ n+1/5 J(η n )dξ = 0, ∀ˆ qˆ∇ q ∈ L2 (Ω)
(3.1)
b Ω
b via the solution of the following pure advection problem ˆ n+2/5 ∈ V (Ω) 2. Find u n n+1 b in Ω × (t , t ) Z Z n b u·v ˆ J(η n )dξ = 0, ˆ·v ˆ J(η )dξ + (ˆ un+1/5 · ∇)ˆ ∂t u b b Ω Ω ˆ (tn ) = u ˆ n+1/5 , u b n+1/5,− × (tn , tn+1 ) , ˆ (t) = u ˆ n+1/5 on Γ u
∀ˆ v ∈ Vb n+1/5,− ,
(3.2) ˆ n+2/5 = u ˆ (tn+1 ). The functional spaces introduced here are and then set u defined as follows
b = V (Ω)
Vb n+1/5,− b n+1/5,− Γ
=
=
b 2 , v|x1 =0 = v|x1 =L } , {v | v ∈ (H 1 (Ω)) b 2, v = 0 {v | v ∈ (H 1 (Ω))
b n+1/5,− , v|x1 =0 = v|x1 =L }, on Γ
b u ˆ n+1/5 · e2 < 0}. {x | x ∈ Γ,
8
Roland Glowinski and Giovanna Guidoboni
JMFM
b and pˆn+3/5 ∈ L20 (Ω) b such that ˆ n+3/5 ∈ V0 (Ω) 3. Find u Z n+3/5 Z ˆ ˆ n+2/5 u −u b un+3/5 ) : D(ˆ b v) J(η n )dξ ˆ J(η n )dξ + µ D(ˆ ̺ ·v ∆t b b Ω Ω Z b ·v b , ˆ J(η n )dξ = 0, ∀ˆ v ∈ V0 (Ω) − pˆn+3/5 ∇ b Ω Z b ·u b . ˆ n+3/5 J(η n )dξ = 0, ∀ˆ qˆ∇ q ∈ L2 (Ω)
(3.3)
b Ω
b via the solution of the following pure advection problem ˆ n+4/5 ∈ V (Ω) 4. Find u n n+1 b in Ω × (t , t )
Z Z n b u·v ˆ J(η n )dξ = 0, ˆ ˆ (wn+3/5 · ∇)ˆ u · v J(η )dξ − ∂ t b b Ω Ω n n+3/5 ˆ ˆ u (t ) = u , b n+3/5,+ × (tn , tn+1 ) , ˆ (t) = u ˆ n+3/5 on Γ u
c n+3/5,+ , ∀ˆ v∈W
(3.4)
ˆ n+4/5 = u ˆ (tn+1 ) and pˆn+4/5 = pˆn+1/5 + pˆn+3/5 . Here and then set u wn+3/5 (ξ1 ) = and
ξ2 n+3/5 n+3/5 u ˆ2 (ξ1 , H) − uˆ1 (ξ1 , H)∂1 η n (ξ1 ) e2 , H
c n+3/5,+ W
=
b n+3/5,+ Γ
=
b 2, v = 0 {v | v ∈ (H 1 (Ω))
(3.5)
b n+3/5,+ , on Γ
v periodic in x1 of period L},
b wn+3/5 · e2 > 0}. {x | x ∈ Γ,
5. Find the new position of the free surface η n+1 ∈ S by solving the following problem in (0, L) × (tn , tn+1 ) Z L Z L n+4/5 n+4/5 (ξ1 , H)) φdξ1 = 0, ∀φ ∈ S, (ξ1 , H)∂1 η − uˆ2 ∂ η φ dξ + (ˆ u1 t 1 0
0
n n η(ξ1 , t ) = η (ξ1 ).
(3.6) b The new physical domain Ω is obtained from the reference domain Ω −1 through the mapping Atn+1 defined as: A−1 tn+1 :
b ∈ R2 Ω
ξ = (ξ1 , ξ2 )
n+1
→ Ωn+1 ∈ R2
x1 = ξ1 → x = (x1 , x2 ) = η n+1 (ξ1 ) x2 = ξ2 . H
(3.7)
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
9
The velocity un+1 and the pressure pn+1 defined on the new domain Ωn+1 are obtained as follows ˆ n+4/5 (x1 , Hx2 /η n+1 (x1 )), un+1 (x1 , x2 ) = u pn+1 (x1 , x2 ) = pˆn+4/5 (x1 , Hx2 /η n+1 (x1 )),
(3.8)
where (x1 , x2 ) ∈ Ωn+1 . b = J−t (η n )∇ and that indeed it is much easier We remark that in (3.1)-(3.4) ∇ to solve the sub-problems (3.1)-(3.3) in their equivalent formulation in the physical domain Ωn (which is fixed from the way the problem was splitted). Moreover, problems (3.1) and (3.3) can be significantly simplified by noticing that Z Z Z µ µ ∇u : ∇v dx + (∇u)t : ∇v dx, (3.9) D(u) : D(u) dx = µ 2 2 Ω(t) Ω(t) Ω(t) and then evaluating the last integral in (3.9) at the previous fractional step. Problems (3.1)-(3.3) are then modified as follows. 1. Find un+1/5 ∈ V0 (Ωn ) and pn+1/5 ∈ L20 (Ωn ) such that Z Z Z un+1/5 − un µ n+1/5 ̺ · vdx + pn+1/5 ∇ · vdx ∇u : ∇vdx − ∆t 2 n n n Ω Ω Ω Z Z Z µ g · vdx − H(η n )n · vdΓn + =s (∇un )t : ∇vdx, ∀v ∈ V0 (Ωn ) , 2 n n n Ω Γ Ω Z n+1/5 2 n q∇ · u dx = 0, ∀q ∈ L (Ω ) . Ωn
(3.10) 2. Find un+2/5 ∈ V (Ωn ) via the solution of the following pure advection problem in Ωn × (tn , tn+1 ) Z Z (un+1/5 · ∇)u · vdx = 0, ∂t u · vdx + Ωn Ωn u(tn ) = un+1/5 , u(t) = un+1/5 on Γn+1/5,− × (tn , tn+1 ) ,
∀v ∈ V n+1/5,− , (3.11)
and then set un+2/5 = u(tn+1 ). The functional spaces introduced here are defined as follows V (Ωn ) = {v | v ∈ (H 1 (Ωn ))2 , v|x1 =0 = v|x1 =L } ,
V n+1/5,− = {v | v ∈ (H 1 (Ωn ))2 , v = 0
on Γn+1/5,− , v|x1 =0 = v|x1 =L },
Γn+1/5,− = {x | x ∈ Γn , un+1/5 · n(x) < 0}.
10
Roland Glowinski and Giovanna Guidoboni
JMFM
3. Find un+3/5 ∈ V0 (Ωn ) and pn+3/5 ∈ L20 (Ωn ) such that Z Z Z un+3/5 − un µ n+3/5 ̺ · vdx + pn+3/5 ∇ · vdx ∇u : ∇vdx − ∆t 2 n n n Ω Ω Ω Z µ =− (∇un+2/5 )t : ∇vdx, ∀v ∈ V0 (Ωn ) , 2 n Ω Z q∇ · un+3/5 dx = 0, ∀q ∈ L2 (Ωn ) . Ωn
(3.12)
ˆ n+3/5
ˆ n+3/5
Then we define u by u is needed to solve problem (3.4).
n+3/5
(ξ1 , ξ2 ) = u
n
ˆ n+3/5 (ξ1 , η (ξ1 )ξ2 /H); u
4. Space discretization At each time step, the subproblems (3.10)-(3.12) are solved in the domain Ωn which is clearly non-polygonal because of the curved portion Γn of the boundary. To approximate velocity and pressure we use here an isoparametric version (discussed in, e.g., [8], Chapter 5) of the Bercovier-Pironneau finite elements method also known as the P1 − iso − P2 and P1 finite elements approximation. This approximation was introduced in [3] and further discussed in, e.g., [18]. To describe the isoparametric finite element approximation of the problem let us consider a bounded domain Ω(t) ∈ R2 with a curved boundary Γ(t), like the one in Figure 1. We introduce a triangulation Th of Ω with discretization step h S and we decompose it as Th = TRh T0h , where TRh
= {T | T ∈ Th , the three edges of T are rectilinear} ,
(4.1)
T0h
= {T | T ∈ Th , T has one curvilinear edge with the two extremities on Γ(t)} .
(4.2) (4.3)
Every rectilinear triangle T ∈ TRh is divided into four sub-triangles KiT , i = 1, 2, 3, 4 by joining the midpoints of its edges. On the other hand, every curved triangle T ∈ T0h is approximated by the quadrilateral Te obtained by joining the midpoints of the two rectilinear edges and the midpoint of the curved one. Then, the quadrilateral Te is divided into four sub-triangles KiT , i = 1, 2, 3, 4 as in the previous case. To construct the velocity spaces we introduce {q | q ∈ C 0 (T ), q|KiT ∈ P1 , ∀i = 1, 2, 3, 4}, if T ∈ TRh Pe2 (T ) = (4.4) {q | q ∈ C 0 (Te), q|KiT ∈ P1 , ∀i = 1, 2, 3, 4}, if T ∈ T0h
where C 0 (T ) (resp. C 0 (Te)) is the space of the functions continuous over T (resp. Te), and P1 is the space of polynomials with degree less or equal to one. Similarly,
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
11
Figure 3. Triangulation of the domain: subdivision of a rectilinear triangle (on the left) and of a curved triangle (on the right). the following three-dimensional subspace of Pe2 (T ) will be useful to construct the pressure spaces:
Pe1 (T ) = {q | q ∈ Pe2 (T ), q(aijT ) = (q(aiT ) + q(ajT ))/2, ∀i, j 1 ≤ i, j ≤ 3, i 6= j}, (4.5) where aiT , i = 1, 2, 3, are the vertices of T , while aijT , with i, j = 1, 2, 3, are the midpoints of the (curved or rectilinear) edges of vertices aiT and ajT . We use the convention aijT = ajiT , ∀i, j, i 6= j. S Thus, we define Te0h = {Te}T ∈T0h and Teh = TRh S Te0h ; the discrete approximation of the domain Ω is given by Ωh = interior of K∈T˜h K and Γh (t) = ∂Ωh (t). Finally we define the pressure and velocity finite element spaces as Ph (Ω(t))
Vh (Ω(t))
=
=
{qh | qh ∈ C 0 (Ωh (t)), qh |T ∈ P1 , ∀T ∈ TRh , qh |Te ∈ Pe1 (T ), ∀T ∈ T0h , qh |x1 =0 = qh |x1 =L } ,
{vh | vh ∈ (C 0 (Ωh (t)))2 , vh |T ∈ (Pe2 (T ))2 , ∀T ∈ TRh , vh |Te ∈ (Pe2 (T ))2 , ∀T ∈ T0h , vh |x1 =0 = vh |x1 =L } .
(4.6)
(4.7)
In order to define the space for the function describing the shape of the free boundary, we need to describe the details of how the triangulations of the domain are obtained. Let us take a uniform structured triangulation of the reference rectanb see the picture on the left in Figure 4. The triangulation Teh of the gular domain Ω, current domain is obtained through the mapping A−1 which actually corresponds t to either a dilation or a compression in the vertical direction proportional to the local height of the fluid layer. Therefore, the discretization in the x1 direction reb is main unchanged for 0 < t < T . Moreover, since the original triangulation in Ω structured and uniform, the number of nodes on the free boundary Γ(t) is equal to the number of nodes on the rigid bottom Σ, and their coordinates differ only in the x2 component. As a consequence, the rigid surface Σ is divided in intervals of equal size. Therefore we can introduce the set Sh = {ηh | ηh ∈ C 0 [0, L], ηh |I ∈ P1 , ∀I ∈ Ih , ηh |x1 =0 = ηh |x1 =L }, where Ih is the set of the edges of the velocity triangulation included in Σ.
(4.8)
12
Roland Glowinski and Giovanna Guidoboni
JMFM
Figure 4. Triangulation of the domain. Triangulation for the pressure in the reference domain (on the left ), triangulation for the pressure in the current domain (on the center ), triangulation for the velocity in the current domain (on the right ).
5. Fully discrete problem Let us introduce the following spaces: V0h (Ωnh ) = {vh | vh ∈ Vh (Ωnh ), vh = 0 on Σh , vh |x1 =0 = vh |x1 =L } Z qh dx = 0, qh |x1 =0 = qh |x1 =L } P0h (Ωnh ) = {qh | qh ∈ Ph (Ωnh ), n+1/5,− Vh n+1/5,−
Γh
c n+3/5,+ W h b n+3/5,+ Γ h
Ωn h
= {vh | vh ∈
Vh (Ωnh ),
n+1/5,−
vh = 0 on Γh
n+1/5
= {x | x ∈ Γnh , uh
· n(x) < 0} n+3/5,+
b h ), vh = 0 on Γ b = {vh | vh ∈ Vh (Ω h n+3/5
bh , w = {x | x ∈ Γ h
, vh |x1 =0 = vh |x1 =L }
, vh |x1 =0 = vh |x1 =L },
· e2 > 0}.
Then, the fully discretized problem reads as follows. n+1/5
n+1/5
1. Find uh ∈ V0h (Ωnh ) and ph ∈ P0h (Ωnh ) such that Z Z Z n+1/5 uh − unh µ n+1/5 n+1/5 ̺ · v dx + ∇u : ∇v dx − ph ∇ · vh dx h h h n n n ∆t 2 Ωh Ωh Ωh Z L Z L n 2 −1/2 ∂1 ηhn (1 + |∂1 ηhn |2 )−1/2 ∂1 v2h dx1 (1 + |∂ η | ) ∂ v dx − s = −s 1 h 1 1h 1 0 0 Z Z µ n t + (∇uh ) : ∇vh dx, ∀vh ∈ V0h (Ωnh ) , g · vh dx − n n 2 Ω Ω h h Z n+1/5 qh ∇ · uh dx = 0, ∀qh ∈ Ph (Ωnh ) , Ωn h
where Remark 1 has been taken into account.
(5.1)
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
13
n+2/5
2. Find uh ∈ Vh (Ωnh ) via the solution of the following pure advection probn lem in Ωh × (tn , tn+1 ) Z Z n+1/5 ∂t uh · vh dx + (uh · ∇)uh · vh dx = 0, n n Ωh Ωh n+1/5 uh (tn ) = uh , n+1/5 n+1/5,− uh (t) = uh on Γh × (tn , tn+1 ) ,
n+1/5,−
∀vh ∈ Vh
, (5.2)
n+2/5
and then set uh = uh (tn+1 ). n+3/5 n+3/5 3. Find uh ∈ V0h (Ωnh ) and ph ∈ P0h (Ωnh ) such that Z Z Z n+3/5 uh − unh µ n+3/5 n+3/5 · vh dx + ∇uh : ∇vh dx − ph ∇ · vh dx ̺ n n n ∆t 2 Ω Ω Ω h h h Z µ n+2/5 t =− (∇uh ) : ∇vh dx, ∀vh ∈ V0h (Ωnh ) , n 2 Ωh Z n+3/5 qh ∇ · uh dx = 0, ∀q ∈ Ph (Ωnh ) . Ωn h
(5.3) b 4. Find ∈ Vh (Ωh ) via the solution of the following pure advection probb h × (tn , tn+1 ) lem in Ω n+4/5 ˆh u
Z Z n+3/5 b ˆh J(η)dξ = 0, ˆ ˆ (wh · ∇)ˆ uh · v u · v J(η)dξ − ∂ h t h b b Ωh Ωh n+3/5 ˆ h (tn ) = u ˆh u , n+3/5 b n+3/5,+ × (tn , tn+1 ) , ˆ h (t) = u ˆh u on Γ h
n+3/5,+
c ∀ˆ vh ∈ W h
,
(5.4) n+4/5 ˆh ˆ h (tn+1 ) and pn+4/5 = pn+1/5 + pn+1/3 . and then set u =u 5. Find the new position of the free surface ηhn+1 ∈ Sh by solving the following problem in (0, L) × (tn , tn+1 )
Z
0
L
∂t ηh φh dξ1 +
Z
n n ηh (ξ1 , t ) = ηh (ξ1 ),
0
L
n+4/5
(ˆ u1h
n+4/5
(ξ1 , H)∂1 ηh − u ˆ2h
(ξ1 , H)) φh dξ1 = 0,
∀φh ∈ Sh ,
(5.5)
14
Roland Glowinski and Giovanna Guidoboni
JMFM
b h trough the The new domain Ωn+1 is obtained from the reference domain Ω h mapping A−1 defined as: tn+1 ,h A−1 tn+1 ,h :
b h ∈ R2 Ω
→ Ωn+1 ∈ R2 h
x1 = ξ1 ξ = (ξ1 , ξ2 ) → x = (x1 , x2 ) = η n+1 (ξ1 ) x2 = h ξ2 . H
(5.6)
The velocity un+1 and the pressure pn+1 defined on the new domain Ωn+1 h h h are obtained as follows n+4/5
ˆh un+1 (x1 , x2 ) = u h
(x1 , ηhn (x1 )x2 /η n+1 (x1 )),
n+4/5
pn+1 (x1 , x2 ) = ph h
(x1 , ηhn (x1 )x2 /ηhn+1 ),
(5.7) (5.8)
where (x1 , x2 ) ∈ Ωn+1 . h The backward Euler’s method has been used to derive the time discretization in (5.1) and (5.3). These problem can be seen as a ”generalized” discrete Stokes problem for which efficient solution methods already exist as shown in, e.g., [8]. The advection step (3.2) has been solved using a wave-like equation method [8, 10, 17]: this approach preserves the hyperbolic nature of advection, it does not introduce numerical dissipation and it is easy to implement. It consists in replacing problem (5.2) by the following wave-like equation one in Ωnh × (tn , tn+1 ): Z Z n+1/5 n+1/5 2 ∂ u · v dx + (uh · ∇)uh · (uh · ∇)vh dx h t h n n Ωh Ωh Z n+1/5 n+1/5,− + (uh · n)(∂t uh · vh )dΓnh = 0, ∀vh ∈ V0h , n+1/5,− Γn \Γ h h n+1/5 uh (tn ) = uh , Z Z n+1/5 n+1/5,− n ∂ u (t ) · v dx = − (un+1/5 · ∇)uh · vh dx, ∀vh ∈ V0h , t h h n n Ωh Ωh n+1/5,− n+1/5 n+1/5,− ∂t uh (tn ) ∈ V0h , uh = uh on Γh × (tn , tn+1 ). (5.9) We have solved problem (5.9) using a second order accurate time discretization scheme which is discussed in, e.g., [8], Chapter 6, and [17]. To solve the trasport problem (5.4) we use again the wave equation approach. ˆ and therefore the We observe that (5.4) does not contain x1 differentiation of u problem reduces to the solution of a family (infinite for the continuous problem, finite for the discrete ones) of transport problems in one space dimension along ˆ is solution of a the vertical direction. Then for ξ1 ∈ [0, L), each component of u
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
transport problem of the following form: ∂ϕ ∂ϕ = 0 on (0, L) × (tn , tn+1 ), − aξ2 ∂t ∂ξ2 ϕ(tn ) = ϕ0 , ϕ(H, t) = b if a > 0, t ∈ (tn , tn+1 ),
15
(5.10)
with a and b constant with respect to ξ2 and t. We observe that any smooth solution of (5.10) verifies: 2 ∂ ϕ ∂ϕ ∂ n n+1 2 ), ∂t2 − a ξ2 ξ2 ξ2 ξ2 = 0 on (0, L) × (t , t ϕ(tn ) = ϕ0 , ∂ϕ ∂ϕ0 , (tn ) = aξ2 (5.11) ∂t ∂ξ2 n n+1 ϕ(H, t) = b if a > 0, t ∈ (t , t ), ∂ϕ ∂ϕ = 0 if a ≤ 0, t ∈ (tn , tn+1 ). aξ2 ∂t − aξ2 ∂ξ 2 ξ2 =H
From a practical point of view, we take advantage of the following variational formulation of (5.10), better suited for finite element treatment. 1. Case a > 0. Z H 2 Z H ∂ ϕ v ∂ϕ ∂v 2 dξ + a dξ2 = 0, ξ2 2 2 ξ ∂t ∂ξ 2 2 ∂ξ2 0 0 ϕ(t0 ) = ϕ0 ,
∀v ∈ D(0, H), a.e. on (t0 , tf ),
∂ϕ0 ∂ϕ , (t0 ) = aξ2 ∂t ∂ξ2 ϕ(H, t) = f,
with D(0, H) = {v | v ∈ C ∞ [0, H], v has compact support in (0, H) }. 2. Case a < 0. Z H 2 Z H ∂ϕ ∂ϕ ∂v ∂ ϕ v 2 dξ2 + a dξ2 − a (H, t)v(H) = 0, ξ2 2 ξ ∂t ∂ξ ∂ξ ∂t 2 2 2 0 0 ∀v ∈ D(0, H), a.e. on (t0 , tf ), ϕ(t0 ) = ϕ0 , ∂ϕ0 ∂ϕ , (t0 ) = aξ2 ∂t ∂ξ2
with V0 = {v | v ∈ C ∞ [0, H], v has compact support in (0, H] }.
(5.12)
(5.13)
16
Roland Glowinski and Giovanna Guidoboni
JMFM
Finally, problem (5.5) is solved using a second order Taylor-Galerkin discretization on (0, L) × (tn , tn+1 ), namely Z L Z L n+4/5 n+4/5 (ηhn+1 − ηhn ) φh dx1 = ∆t (u2h − u1h ∂1 ηhn ) φh dx1 0 0 Z (∆t)2 L n+4/5 n+4/5 (5.14) u1h ∂1 u2h φh dx1 − 2 0 Z (∆t)2 L n+4/5 n+4/5 n+4/5 − u1h ∂1 ηhn (∂1 u1h φh + u1h ∂1 φh ) dx1 , 2 0 where the velocity is computed at Γnh . Remark 2. For a general approach concerning the approximation of free surface problems for incompressible viscous fluids with surface tension see, e.g., [4].
6. Numerical experiments Problem (2.4)-(2.10) always admits a steady solution. Indeed, assuming that u = U (x2 )e1 , with e1 = (1, 0), p = P (x2 ) and η(x1 , t) ≡ H, we get U=
̺g sin β (2Hx2 − x22 ), 2µ
P = −̺g cos β(x2 − H) in (0, L) × (0, H). (6.1)
Nishida et al. in [15] proved the existence of time periodic solutions bifurcating from the basic flow (6.1) for Reynolds numbers large enough. The Reynolds number Re is defined as Re = HUmax ̺/µ, where H is the average height of the fluid layer and Umax is the velocity of the basic flow at the free surface, obtained from (6.1) with x2 = H, namely H 2 ̺g sin β Umax = , (6.2) 2µ implying that H 3 ̺2 g sin β Re = . (6.3) 2µ2 Our goal is to explore the features of the bifurcating solutions using the numerical algorithm described in the previous sections. For the time discretization we have taken ∆t = 2.5 · 10−4 . As mentioned in Section 5, we solve the advection steps (5.2) and (5.4) using a wave-like equation method which requires the specification of time substeps: here we have chosen both of them equal to ∆t/5. In our numerical experiments, the dimensions of the reference rectangular domain are H = 0.15 and L = 3, the density of the fluid is ̺ = 1, and the angle of inclination is β = 45◦ . By direct simulation, we estimated the critical Reynolds number Rec as a function of the surface tension. In particular, we considered s = 1, 2, 3. The results are reported in Table 1. It is clear that Rec is an increasing function of s, which means that the surface tension has a stabilizing effect, as shown in [13].
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
17
range of Rec range for µc surface tension s s = 1 1.68 < Rec < 1.98 1/12 < µc < 1/13 s = 2 2.29 < Rec < 2.63 1/14 < µc < 1/15 3.00 < Rec < 3.38 1/16 < µc < 1/17 s = 3 Table 1. Range of the critical Reynolds number Rec and the correspondent viscosities for several values of the surface tension s.
The system reaches a steady state whenever the Reynolds number is below the critical one. In the simulations, the solution is considered a steady state when the following condition is satisfied: kun − un−1 kL2 (Ωn ) < 10−7 . kun kL2 (Ωn )
(6.4)
In Figure 5 we compare the analytical velocity profile given by (6.1) and the one obtained numerically with our algorithm (here s = 1, µ = 1/11, and Re = 1.42). The agreement between the two profiles is indeed very satisfactory. When the Reynolds number exceeds the critical threshold, periodic solutions establish. In Figure 8 this phenomenon is tracked by the vertical oscillations of the middle point on the free surface (here s = 1, µ = 1/15, and Re = 2.63). In Figure 9 we have plotted five pictures of the streamlines contours showing the evolution of the solution at Re = 2.63 and s = 1 during one whole period. The wavelength of the waves on the free surface is very long at the early stages of the onset of the instability, as shown in Figure 6. As time goes by, the system enters a periodic regime and the surface waves develop into a profile whose shape, amplitude and frequency depend on the Reynolds number and on the surface tension. Numerical results for s = 1 and different Reynolds numbers are reported in Table 2, and they show that both the amplitude and the frequency of the waves are increasing functions of Re. Moreover, the shape of the waves steepen as the Reynolds number increases, as shown in Figure 7. We repeated the computations for s = 1 and Re = 1.98 using different mesh sizes and time steps. The results are reported in Table 3 and they show a very good agreement in terms of amplitude, frequency and period of the periodic oscillations. Animations of the numerical simulations described in this article are available at http://www.math.uh.edu/∼gio/
References [1] M. Behr, F. Abraham, Free-Surface Flow Simulations in the Presence of Inclined Walls, Comput. Methods Appl. Mech. and Engrg. 191 (2002) 5467-5483 [2] T.B. Benjamin Wave formation in laminar flow down an inclined plane, J. Fluid Mech. 2 (1957) 554-574
18
Roland Glowinski and Giovanna Guidoboni
JMFM
0.14
0.12
x2
0.1
0.08
0.06
0.04
0.02
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
Figure 5. Comparison between the velocity profile at the steady state computed analytically (solid line) and numerically (×). Here s = 1, µ = 1/11, and Re = 1.42. Re amplitude period frequency 1.98 0.03215 1.63150 3.85117 2.63 0.06956 1.41775 4.43180 0.13927 1.15650 5.43293 3.79 Table 2. Amplitude and frequency for the oscillations of the middle point on the free surface for s = 1 and several Reynolds number beyond the Hopf bifurcation. amplitude period frequency (hp , ∆t) (1/100, 2.5 × 10−4 ) 0.03215 1.63150 3.85117 (1/64, 2.5 × 10−4 ) 0.03210 1.63275 3.84822 0.03176 1.6291 3.85691 (1/100, 1.25 × 10−4 ) Table 3. Amplitude and frequency for the oscillations of the middle point on the free surface for s = 1, µ = 1/13 and Re = 1.98 for different combinations of mesh sizes and time steps.
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
19
0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0
0
0.5
1
1.5
2
2.5
3
Figure 6. The instability begins with waves characterized by very long wavelengths. Here s = 1, µ = 1/15, and Re = 2.63. [3] M. Bercovier, O. Pironneau, Error estimates for finite element method solution of the Stokes problem in primitive variables, Numer. Math. 33 (1979) 211-224 [4] A. Caboussat, A numerical method for the simulation of free surface flows with surface tension, Computers and Fluids 35-10 (2006) 1205-1216 [5] R.W. Chin, F.H. Abernathy, R. Bertschy, Gravity and shear wave stability of free surface flows. Part1. Numerical calculations, J. Fluid Mech. 168 (1986) 501-513 [6] W. Dettmer, D. Peri´c, A computational framework for free surface accounting for surface tension, Comput. Methods Appl. Mech. Engrg. 195 (2006) 3038-3071 [7] J. Donea, Arbitrary Lagrangian Eulerian Methods, in Computational Methods for Transient Analysis, Vol. I of Computational Methods in Mechanics. North-Holland, Elsevier (1983) [8] R. Glowinski, Finite Element Methods for Incompressible Viscous Flow, in Handbook of Numerical Analysis, Vol. IX, P.G. Ciarlet and J.-L. Lions eds., North-Holland, Amsterdam (2003) [9] R. Glowinski, L.H. Ju´ arez, Finite element method and operator splitting for a timedependent viscous incompressible free surface flow, Comput. Fluid Dyn. J. 12-3 (2003) 459-468 [10] R. Glowinski, G. Guidoboni, T.W. Pan, Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity, J. Comput. Phys. 216-1 (2006) 76-91
20
Roland Glowinski and Giovanna Guidoboni
JMFM
0.24
0.22
0.2
0.18
0.16
0.14
0.12
0.1
0.08
0
0.5
1
1.5
2
2.5
3
Figure 7. Shape of the waves on the free surface in the periodic regime for Re = 1.98 and µ = 1/13 (solid line), Re = 2.63 and µ = 1/15 (dashed line), Re = 3.79 and µ = 1/18 (dashed-dotted line), when s = 1. [11] D.Rh. Gwynllyw, D.H. Peregrine, Numerical Simulation of Stokes flow down an inclined plane, Proc. R. Soc. Lond. A 452-1946 (1996) 543-565 [12] T.J. Hughes, W.K. Liu, T.K. Zimmermann, Lagrangian-eulerian finite element formulation for incompressible viscous flows, Comp. Meth. Appl. Mech. Engrg. 29 (1981) 329-349 [13] J. Liu, J.D. Paul, J.P. Gollub, Measurements of the primary instabilites of film flows, J. Fluid Mech. 250 (1993) 69-101 [14] G.I. Marchuk, Splitting and alternate direction methods, in Handbook of Numerical Analysis, Vol. I, P.G. Ciarlet and J.-L. Lions eds., North-Holland, Amsterdam (1990) [15] T. Nishida, Y. Teramoto, H. Yoshihara, Hopf bifurcation in viscous incompressible flow down an inclined plane, J. Math. Fluid Mech. 7 (2005) 29-71 [16] M. Padula, V.A. Solonnikov, On Rayleigh-Taylor stability, Ann. Univ. Ferrara Sez. VII (N.S.) 46 (2000), 307-336 [17] T.W. Pan and R. Glowinski, A projection/wave-like equation method for the numerical simulation of incompressible viscous fluid flow modeled by the Navier-Stokes equations , Comput. Fluid Dyn. J. 9-2 (2000) 28-42 [18] O. Pironneau, Finite Element Methods for Fluids, J. Wiley, Chichester (1989)
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
21
0.18
0.17
0.16
0.15
0.14
0.13
0.12 0
5
10
15
20
25
30
35
Figure 8. Variation in time of height of the middle point of the free surface (s = 1, µ = 1/15 and Re = 2.63).
Figure 9. Streamlines of the solution during one period, obtained for s = 1, µ = 1/15 and Re = 2.63. Acknowledgment The authors would like to thank Professors H. Juarez, M. Padula, T.W. Pan, M. Picasso and J. Rappaz for helpful comments and suggestions. The support of
22
Roland Glowinski and Giovanna Guidoboni
JMFM
NSF (Grant DMS-9973318, 0209066, and CCR-9902035) and DOE/LACSI (Grant R71700K-292-000-99) is also acknowledged. Roland Glowinski Department of Mathematics, University of Houston, Houston, Texas 77204-3476, USA and Laboratoire Jacques-Louis Lions, Universit´e P. et M. Curie, 4 Place Jussieu, 75005, Paris, France e-mail: [email protected] Giovanna Guidoboni Department of Mathematics, University of Houston, Houston, Texas 77204-3476, USA e-mail: [email protected]
auser Verlag Basel/Switzerland
Journal of Mathematical Fluid Mechanics
Hopf bifurcation in viscous incompressible flow down an inclined plane: a numerical approach Roland Glowinski and Giovanna Guidoboni Abstract. In this article we investigate numerically a Hopf bifurcation phenomenon for a viscous incompressible flow down an inclined plane. This problem has been discussed by Nishida et al. who proved the existence of periodic solutions bifurcating from the steady flow. Using a computational methodology based on finite elements for the space discretization and on operator splitting for the time discretization, we have been able to reproduce the results predicted by Nishida et al. Mathematics Subject Classification (2000). 70K50, 76B45, 76M10. Keywords. free surface, surface tension, Hopf bifurcation, operator splitting, finite elements.
1. Introduction In [15], Nishida, Teramoto and Yoshihara proved the existence of a Hopf bifurcation for the flow of an incompressible viscous fluid down an inclined plane. The main goal of this article is to investigate numerically the same flow problem in order to precise at which value of the Reynolds number the Hopf bifurcation is occurring, and to visualize the shape of the past-bifurcation waves. The computational methodology discussed in this article relies on the following components: b (i) As in [9] we map the two dimensional flow region Ω(t) on a rectangular one Ω (actually the main result in [15] is proved using a similar transformation, which is quite classical in this context). b and we use a time discretization by (ii) We formulate the flow problem in Ω operator splitting to reduce the solution of the original flow problem to that of simpler sub-problems. Actually, it is easier to solve some of the sub-problems in
2
Roland Glowinski and Giovanna Guidoboni
JMFM
the physical flow region, and we take advantage of this simplification each time it is possible. (iii) We combine the time-splitting scheme with an isoparametric Bercoveir-Pironneau finite element approximation of the Navier-Stokes equations (see e.g. [3], [8], [9], [10], [18]) to fully discretize the free surface problem. (iv) The wave-like equation methodology discussed in, e.g., [8], [9], [10], and [17], is used to handle the pure advection problems resulting from the time-splitting and from the domain transformation. The resulting methodology is modular, relatively easy to implement, and it introduces very little numerical dissipation, an important property if one wishes to identify accurately the critical Reynolds number at which bifurcation occurs. We took advantage of this property in, e.g., [10] to study other Hopf bifurcation phenomena. The computed results are in very good agreement with the experimental data reported in [13]. The flow problem in this article is a basic one. It is not surprising therefore that it has attracted the attention of many investigators. Let us mention, among many others, [1], [2], [5], [6], [11] and [16]. The content of the remainder of this article is as follows. In Section 2 we describe the mathematical formulation of the problem. In Section 3 we discuss the time discretization by operator splitting. In Sections 4 and 5 the splitting scheme of Section 3 is combined with the Bercovier-Pironneau finite element approximation of the Navier-Stokes equations to derive the full discretization of the problem. The results of the numerical experiments are presented in Section 6.
2. Mathematical formulation An incompressible viscous fluid flows down an inclined plane under the action of gravity. Let us consider a Cartesian frame of reference with the x1 -axis (resp. x2 -axis) along (resp. perpendicular to) the plane which is inclined at an angle β, as shown in Figure 1. We assume that the free surface Γ(t) admits a cartesian representation at each instant of time so that it can be described by the timedependent function η(x1 , t). Therefore, the domain occupied by the fluid is Ω(t) = {x = (x1 , x2 ) ∈ R2 : x1 ∈ (0, L), x2 ∈ (0, η(x1 , t))},
(2.1)
the rigid lower boundary is Σ = {x = (x1 , x2 ) ∈ R2 : x1 ∈ (0, L), x2 = 0},
(2.2)
and the upper free surface is Γ(t) = {x = (x1 , x2 ) ∈ R2 : x1 ∈ (0, L), x2 = η(x1 , t)}.
(2.3)
The fluid flow is modeled by the Navier-Stokes equations ∇·u =0 ̺[∂t u + (u · ∇)u] = ∇ · σ + ̺g,
(2.4)
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
x
3
2
Γ(t)
η(x1,t)
0
β L
x
1
Figure 1. Sketch of the geometry of the problem. where u = (u1 , u2 ) is the velocity of the fluid, ̺ is the fluid density, σ is the stress tensor defined as σ = −pI + 2µD(u), where p is the pressure, µ is the viscosity and D(u) = (∇u + (∇u)t )/2 is the rate-of-strain tensor, where the superscript t means transpose. The external force is due to the gravity and it is given by ̺g = ̺g(sin β, − cos β), where g is the acceleration of gravity. Here ∂t denotes the partial derivative with respect to time, while ∇ = (∂1 , ∂2 ) is the gradient with respect to the spatial variables. Equations (2.4) hold in the time-dependent domain Ω(t) defined in (2.1), and for 0 < t < T . No-slip condition is assumed for the velocity at the rigid bottom, and therefore we have u=0
on Σ.
(2.5)
The free boundary Γ(t) is a material surface, therefore we prescribe the following kinematic condition ∂t η + u1 ∂1 η − u2 = 0 on Γ(t), (2.6) to ensure that the fluid particles do not cross the surface. An additional condition is needed to establish the balance of stress σn = sH(η(x1 , t))n on Γ(t), (2.7) p 2 where n = (−∂1 η, 1)/ 1 + |∂1 η| is the outward unit normal vector, s is the surface tension (which is assumed to be positive and constant), and H(η(x1 , t)) = ∂12 η/(1 + |∂1 η|2 )3/2 is the curvature. We assume periodicity in the x1 -direction of given period L. The incompressibility of the fluid and the periodic boundary conditions imply that the area of the fluid domain must be conserved. Indeed we have: Z Z L Z Z d L η dx1 . (2.8) ∂t η dx1 = u · n dΓ(t) = ∇ · ud x = 0= dt 0 0 Γ(t) Ω(t) This means that
Z
0
L
η(x1 , t) dx1 = constant = HL,
(2.9)
4
Roland Glowinski and Giovanna Guidoboni
JMFM
where H is the average height of the fluid layer. To complete the formulation of the problem, initial conditions should be prescribed for the velocity of the fluid and the shape of the domain: u(x, 0) = u0 (x)
with ∇ · u0 = 0, Z L η(x1 , 0) = η0 (x1 ) with η0 (x1 ) dx1 = HL.
(2.10)
0
In order to write the variational formulation of the problem (2.4)-(2.10) we introduce the following functional spaces: V0 (Ω(t))
=
L20 (Ω(t))
=
S
{v | v ∈ (H 1 (Ω(t)))2 , v = 0 on Σ, v|x1 =0 = v|x1 =L } , (2.11) Z qdx = 0} , (2.12) {q | q ∈ L2 (Ω(t)), Ω(t)
=
1
{η | η ∈ H (0, L), η|x1 =0 = η|x1 =L } .
(2.13)
Now, we multiply equations (2.4)1 and (2.4)2 by the test functions q ∈ L2 (Ω(t)) and v ∈ V0 (Ω(t)), respectively, and we use the notation ϕ(t) to denote the function x → ϕ(x, t). The resulting equations are integrated over Ω(t) to obtain the variational formulation of (2.4)-(2.10): For a.e. t > 0, find u(t) ∈ V0 (Ω(t)), p(t) ∈ L20 (Ω(t)), η(t) ∈ S such that Z Z Z ̺ ∂t u · v dx + ̺ (u · ∇)u · v dx + 2µ D(u) : D(v) dx Ω(t) Ω(t) Ω(t) Z Z − p∇ · v dx = ̺ g · v dx Ω(t) Ω(t) Z +s H(η(x1 , t))n · v dΓ(t), ∀v ∈ V0 (Ω(t)), Γ(t) Z q∇ · u dx = 0, ∀q ∈ L2 (Ω(t)), Ω(t) Z L Z L ∂t η φ dx1 + (u1 |Γ(t) ∂1 η − u2 |Γ(t) ) φ dx1 = 0, ∀φ ∈ S, 0
(2.14)
0
u(x, 0) = u0 (x), with ∇ · u0 = 0, Z L η(x1 , 0) = η0 (x1 ) with η0 (x1 ) dx1 = HL. 0
Remark 1. There is no contradiction between the choice of S given by (2.13) and the fact that H(η) includes ∂12 η. Indeed, by integration by parts we get Z Z L s H(η)v1 n1 dΓ(t) = −s (1 + |∂1 η|2 )−1/2 ∂1 v1 dx1 , (2.15) Γ(t)
0
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
5
At
H
0
L
0
Ω (t)
L
Ω
−1
At
Figure 2. At maps the the current domain Ω(t) into the referb ence domain Ω. s
Z
H(η)v2 n2 dΓ(t) = −s
Z
L
∂1 η(1 + |∂1 η|2 )−1/2 ∂1 v2 dx1 ,
(2.16)
0
Γ(t)
where the x1 -differentiation on η is transferred on the two components of the test function v. Problem (2.14) is defined on the domain Ω(t) which changes in time and it is not known a priori. Therefore we are going to reformulate problem (2.14) into a fixed reference domain. Let At be a family of mappings which at each time t ∈ (0, T ) maps the b = (0, L) × (0, H) defined by current domain Ω(t) into the reference domain Ω At :
Ω(t) ∈ R2
b ∈ R2 → Ω
x = (x1 , x2 ) → ξ = (ξ1 , ξ2 ) = At (x) =
ξ1 = x1 ξ2 =
H x2 , η(x1 , t)
(2.17)
as sketched in Figure 2. We observe that the free surface Γ(t) is mapped into ˆ = {ξ ∈ R2 | ξ1 ∈ (0, L), ξ2 = H}. Γ The Jacobian matrix J of this transformation 1 0 J = J(η) = ξ2 H − ∂1 η(ξ1 , t) η(ξ1 , t) η(ξ1 , t)
(2.18)
and its determinant J are given by ,
J = J(η) =
H . η(ξ1 , t)
(2.19)
It is clear that this transformation is well defined as long as η(ξ1 , t) > 0 or, in other words, as long as the free surface does not touch the bottom. Let f = f (x, t) be a function defined on Ω(t) and fˆ = fˆ(ξ, t) = f (A−1 t (ξ), t) the correspondent
6
Roland Glowinski and Giovanna Guidoboni
b Due to the change of variables, we have function defined on Ω. Z Z fˆ(ξ, t) Jdξ. f (x, t) dx = b Ω
Ω(t)
Moreover, it follows from the chain rule that b fˆ, ∂t f = ∂t fˆ − w · ∇
JMFM
(2.20)
(2.21)
b = J−t ∇ and the domain velocity w is given by where ∇ ξ2 w(ξ1 , t) = ∂t x = ∂t η(ξ1 , t)e2 H (2.22) ξ2 = u2 (ξ1 , η(ξ1 , t), t) − u1 (ξ1 , η(ξ1 , t), t)∂1 η(ξ1 , t) e2 . H The last step is a consequence of the kinematic condition (2.6) and we used the notation e2 = (0, 1). Using (2.20) and (2.21) , problem (2.14) becomes: b pˆ(t) ∈ L2 (Ω), b η(t) ∈ S, such that ˆ (t) ∈ V0 (Ω), For a.e. t > 0, find u 0
̺
Z b u·v b u·v ˆ J(η)dξ ˆ J(η)dξ + ̺ (ˆ ˆ·v ˆ J(η)dξ − ̺ (w · ∇)ˆ u · ∇)ˆ ∂t u b b b Ω Ω Ω Z Z Z ˆ ·v b u) : D(ˆ b v) J(η)dξ − ˆ·v ˆ J(η)dξ ˆ J(η)dξ = ̺ g pˆ∇ +2µ D(ˆ Z
Z
b Ω
Z
b Ω
b Ω
b H(η) vˆ2 (ξ1 , H, t) − vˆ1 (ξ1 , H, t)∂1 η dξ1 , ∀ˆ v ∈ V0 (Ω), +s 0 Z b ·u b ˆ J(η)dξ = 0, ∀ˆ qˆ∇ q ∈ L2 (Ω), b Ω Z L Z L (ˆ u1 (ξ1 , H, t)∂1 η − uˆ2 (ξ1 , H, t)) φ dξ1 = 0, ∀φ ∈ S, ∂t η φ dξ1 + 0
L
(2.23)
0
b ·u ˆ (ξ, 0) = u ˆ 0 (ξ), with ∇ ˆ 0 J(η0 ) = 0, u Z L η(ξ1 , 0) = η0 (ξ1 ) with η0 (ξ1 ) dξ1 = HL. 0
b u) = (∇ˆ b u + (∇ˆ b u)t )/2. Here we used the notation D(ˆ
3. Time discretization The numerical solution of problem (2.23) contains three main difficulties: (a) the incompressibility condition and the related unknown pressure, (b) an advection term associated to the velocity of the fluid, (c) an advection term associated to the velocity of the domain. The difficulty associated with the unknown location of the free boundary has been removed through the transformation defined by (2.17). A time discretization by operator splitting allows to decouple the difficulties of the problem and to treat each of them with appropriate algorithms. Moreover it
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
7
makes possible to update the location of the free boundary as a sub-step, avoiding costly sub-iterations between the location of the free surface and the solution of the Navier-Stokes equations. Among the many possible versions, we advocate a very simple operator splitting method which is first-order accurate [14]. The low order accuracy is compensated, actually, by easy implementation and less cost in computational time. ˆn = u ˆ (tn ), pˆn = pˆ(tn ), Let ∆t be the time discretization step and tn = n∆t, u n n n ˆ and η being known, the scheme consists of solving the η = η(t ). For n ≥ 0, u following problems. n
b and pˆn+1/5 ∈ L2 (Ω) b such that ˆ n+1/5 ∈ V0 (Ω) 1. Find u 0
Z Z n+1/5 ˆ ˆn u −u n b un+1/5 ) : D(ˆ b v)J(η n )dξ ˆ · vJ(η )dξ + µ D(ˆ ̺ ∆t b b Ω Ω Z Z n+1/5 b n ˆ·v ˆ J(η n )dξ ˆ g p ˆ ∇ · v J(η )dξ = ̺ − b b Ω Ω Z L n n b , +s dξ1 , ∀ˆ v ∈ V0 (Ω) H(η ) v ˆ (ξ , H, t) − v ˆ (ξ , H, t)∂ η 2 1 1 1 1 0 Z b ·u b . ˆ n+1/5 J(η n )dξ = 0, ∀ˆ qˆ∇ q ∈ L2 (Ω)
(3.1)
b Ω
b via the solution of the following pure advection problem ˆ n+2/5 ∈ V (Ω) 2. Find u n n+1 b in Ω × (t , t ) Z Z n b u·v ˆ J(η n )dξ = 0, ˆ·v ˆ J(η )dξ + (ˆ un+1/5 · ∇)ˆ ∂t u b b Ω Ω ˆ (tn ) = u ˆ n+1/5 , u b n+1/5,− × (tn , tn+1 ) , ˆ (t) = u ˆ n+1/5 on Γ u
∀ˆ v ∈ Vb n+1/5,− ,
(3.2) ˆ n+2/5 = u ˆ (tn+1 ). The functional spaces introduced here are and then set u defined as follows
b = V (Ω)
Vb n+1/5,− b n+1/5,− Γ
=
=
b 2 , v|x1 =0 = v|x1 =L } , {v | v ∈ (H 1 (Ω)) b 2, v = 0 {v | v ∈ (H 1 (Ω))
b n+1/5,− , v|x1 =0 = v|x1 =L }, on Γ
b u ˆ n+1/5 · e2 < 0}. {x | x ∈ Γ,
8
Roland Glowinski and Giovanna Guidoboni
JMFM
b and pˆn+3/5 ∈ L20 (Ω) b such that ˆ n+3/5 ∈ V0 (Ω) 3. Find u Z n+3/5 Z ˆ ˆ n+2/5 u −u b un+3/5 ) : D(ˆ b v) J(η n )dξ ˆ J(η n )dξ + µ D(ˆ ̺ ·v ∆t b b Ω Ω Z b ·v b , ˆ J(η n )dξ = 0, ∀ˆ v ∈ V0 (Ω) − pˆn+3/5 ∇ b Ω Z b ·u b . ˆ n+3/5 J(η n )dξ = 0, ∀ˆ qˆ∇ q ∈ L2 (Ω)
(3.3)
b Ω
b via the solution of the following pure advection problem ˆ n+4/5 ∈ V (Ω) 4. Find u n n+1 b in Ω × (t , t )
Z Z n b u·v ˆ J(η n )dξ = 0, ˆ ˆ (wn+3/5 · ∇)ˆ u · v J(η )dξ − ∂ t b b Ω Ω n n+3/5 ˆ ˆ u (t ) = u , b n+3/5,+ × (tn , tn+1 ) , ˆ (t) = u ˆ n+3/5 on Γ u
c n+3/5,+ , ∀ˆ v∈W
(3.4)
ˆ n+4/5 = u ˆ (tn+1 ) and pˆn+4/5 = pˆn+1/5 + pˆn+3/5 . Here and then set u wn+3/5 (ξ1 ) = and
ξ2 n+3/5 n+3/5 u ˆ2 (ξ1 , H) − uˆ1 (ξ1 , H)∂1 η n (ξ1 ) e2 , H
c n+3/5,+ W
=
b n+3/5,+ Γ
=
b 2, v = 0 {v | v ∈ (H 1 (Ω))
(3.5)
b n+3/5,+ , on Γ
v periodic in x1 of period L},
b wn+3/5 · e2 > 0}. {x | x ∈ Γ,
5. Find the new position of the free surface η n+1 ∈ S by solving the following problem in (0, L) × (tn , tn+1 ) Z L Z L n+4/5 n+4/5 (ξ1 , H)) φdξ1 = 0, ∀φ ∈ S, (ξ1 , H)∂1 η − uˆ2 ∂ η φ dξ + (ˆ u1 t 1 0
0
n n η(ξ1 , t ) = η (ξ1 ).
(3.6) b The new physical domain Ω is obtained from the reference domain Ω −1 through the mapping Atn+1 defined as: A−1 tn+1 :
b ∈ R2 Ω
ξ = (ξ1 , ξ2 )
n+1
→ Ωn+1 ∈ R2
x1 = ξ1 → x = (x1 , x2 ) = η n+1 (ξ1 ) x2 = ξ2 . H
(3.7)
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
9
The velocity un+1 and the pressure pn+1 defined on the new domain Ωn+1 are obtained as follows ˆ n+4/5 (x1 , Hx2 /η n+1 (x1 )), un+1 (x1 , x2 ) = u pn+1 (x1 , x2 ) = pˆn+4/5 (x1 , Hx2 /η n+1 (x1 )),
(3.8)
where (x1 , x2 ) ∈ Ωn+1 . b = J−t (η n )∇ and that indeed it is much easier We remark that in (3.1)-(3.4) ∇ to solve the sub-problems (3.1)-(3.3) in their equivalent formulation in the physical domain Ωn (which is fixed from the way the problem was splitted). Moreover, problems (3.1) and (3.3) can be significantly simplified by noticing that Z Z Z µ µ ∇u : ∇v dx + (∇u)t : ∇v dx, (3.9) D(u) : D(u) dx = µ 2 2 Ω(t) Ω(t) Ω(t) and then evaluating the last integral in (3.9) at the previous fractional step. Problems (3.1)-(3.3) are then modified as follows. 1. Find un+1/5 ∈ V0 (Ωn ) and pn+1/5 ∈ L20 (Ωn ) such that Z Z Z un+1/5 − un µ n+1/5 ̺ · vdx + pn+1/5 ∇ · vdx ∇u : ∇vdx − ∆t 2 n n n Ω Ω Ω Z Z Z µ g · vdx − H(η n )n · vdΓn + =s (∇un )t : ∇vdx, ∀v ∈ V0 (Ωn ) , 2 n n n Ω Γ Ω Z n+1/5 2 n q∇ · u dx = 0, ∀q ∈ L (Ω ) . Ωn
(3.10) 2. Find un+2/5 ∈ V (Ωn ) via the solution of the following pure advection problem in Ωn × (tn , tn+1 ) Z Z (un+1/5 · ∇)u · vdx = 0, ∂t u · vdx + Ωn Ωn u(tn ) = un+1/5 , u(t) = un+1/5 on Γn+1/5,− × (tn , tn+1 ) ,
∀v ∈ V n+1/5,− , (3.11)
and then set un+2/5 = u(tn+1 ). The functional spaces introduced here are defined as follows V (Ωn ) = {v | v ∈ (H 1 (Ωn ))2 , v|x1 =0 = v|x1 =L } ,
V n+1/5,− = {v | v ∈ (H 1 (Ωn ))2 , v = 0
on Γn+1/5,− , v|x1 =0 = v|x1 =L },
Γn+1/5,− = {x | x ∈ Γn , un+1/5 · n(x) < 0}.
10
Roland Glowinski and Giovanna Guidoboni
JMFM
3. Find un+3/5 ∈ V0 (Ωn ) and pn+3/5 ∈ L20 (Ωn ) such that Z Z Z un+3/5 − un µ n+3/5 ̺ · vdx + pn+3/5 ∇ · vdx ∇u : ∇vdx − ∆t 2 n n n Ω Ω Ω Z µ =− (∇un+2/5 )t : ∇vdx, ∀v ∈ V0 (Ωn ) , 2 n Ω Z q∇ · un+3/5 dx = 0, ∀q ∈ L2 (Ωn ) . Ωn
(3.12)
ˆ n+3/5
ˆ n+3/5
Then we define u by u is needed to solve problem (3.4).
n+3/5
(ξ1 , ξ2 ) = u
n
ˆ n+3/5 (ξ1 , η (ξ1 )ξ2 /H); u
4. Space discretization At each time step, the subproblems (3.10)-(3.12) are solved in the domain Ωn which is clearly non-polygonal because of the curved portion Γn of the boundary. To approximate velocity and pressure we use here an isoparametric version (discussed in, e.g., [8], Chapter 5) of the Bercovier-Pironneau finite elements method also known as the P1 − iso − P2 and P1 finite elements approximation. This approximation was introduced in [3] and further discussed in, e.g., [18]. To describe the isoparametric finite element approximation of the problem let us consider a bounded domain Ω(t) ∈ R2 with a curved boundary Γ(t), like the one in Figure 1. We introduce a triangulation Th of Ω with discretization step h S and we decompose it as Th = TRh T0h , where TRh
= {T | T ∈ Th , the three edges of T are rectilinear} ,
(4.1)
T0h
= {T | T ∈ Th , T has one curvilinear edge with the two extremities on Γ(t)} .
(4.2) (4.3)
Every rectilinear triangle T ∈ TRh is divided into four sub-triangles KiT , i = 1, 2, 3, 4 by joining the midpoints of its edges. On the other hand, every curved triangle T ∈ T0h is approximated by the quadrilateral Te obtained by joining the midpoints of the two rectilinear edges and the midpoint of the curved one. Then, the quadrilateral Te is divided into four sub-triangles KiT , i = 1, 2, 3, 4 as in the previous case. To construct the velocity spaces we introduce {q | q ∈ C 0 (T ), q|KiT ∈ P1 , ∀i = 1, 2, 3, 4}, if T ∈ TRh Pe2 (T ) = (4.4) {q | q ∈ C 0 (Te), q|KiT ∈ P1 , ∀i = 1, 2, 3, 4}, if T ∈ T0h
where C 0 (T ) (resp. C 0 (Te)) is the space of the functions continuous over T (resp. Te), and P1 is the space of polynomials with degree less or equal to one. Similarly,
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
11
Figure 3. Triangulation of the domain: subdivision of a rectilinear triangle (on the left) and of a curved triangle (on the right). the following three-dimensional subspace of Pe2 (T ) will be useful to construct the pressure spaces:
Pe1 (T ) = {q | q ∈ Pe2 (T ), q(aijT ) = (q(aiT ) + q(ajT ))/2, ∀i, j 1 ≤ i, j ≤ 3, i 6= j}, (4.5) where aiT , i = 1, 2, 3, are the vertices of T , while aijT , with i, j = 1, 2, 3, are the midpoints of the (curved or rectilinear) edges of vertices aiT and ajT . We use the convention aijT = ajiT , ∀i, j, i 6= j. S Thus, we define Te0h = {Te}T ∈T0h and Teh = TRh S Te0h ; the discrete approximation of the domain Ω is given by Ωh = interior of K∈T˜h K and Γh (t) = ∂Ωh (t). Finally we define the pressure and velocity finite element spaces as Ph (Ω(t))
Vh (Ω(t))
=
=
{qh | qh ∈ C 0 (Ωh (t)), qh |T ∈ P1 , ∀T ∈ TRh , qh |Te ∈ Pe1 (T ), ∀T ∈ T0h , qh |x1 =0 = qh |x1 =L } ,
{vh | vh ∈ (C 0 (Ωh (t)))2 , vh |T ∈ (Pe2 (T ))2 , ∀T ∈ TRh , vh |Te ∈ (Pe2 (T ))2 , ∀T ∈ T0h , vh |x1 =0 = vh |x1 =L } .
(4.6)
(4.7)
In order to define the space for the function describing the shape of the free boundary, we need to describe the details of how the triangulations of the domain are obtained. Let us take a uniform structured triangulation of the reference rectanb see the picture on the left in Figure 4. The triangulation Teh of the gular domain Ω, current domain is obtained through the mapping A−1 which actually corresponds t to either a dilation or a compression in the vertical direction proportional to the local height of the fluid layer. Therefore, the discretization in the x1 direction reb is main unchanged for 0 < t < T . Moreover, since the original triangulation in Ω structured and uniform, the number of nodes on the free boundary Γ(t) is equal to the number of nodes on the rigid bottom Σ, and their coordinates differ only in the x2 component. As a consequence, the rigid surface Σ is divided in intervals of equal size. Therefore we can introduce the set Sh = {ηh | ηh ∈ C 0 [0, L], ηh |I ∈ P1 , ∀I ∈ Ih , ηh |x1 =0 = ηh |x1 =L }, where Ih is the set of the edges of the velocity triangulation included in Σ.
(4.8)
12
Roland Glowinski and Giovanna Guidoboni
JMFM
Figure 4. Triangulation of the domain. Triangulation for the pressure in the reference domain (on the left ), triangulation for the pressure in the current domain (on the center ), triangulation for the velocity in the current domain (on the right ).
5. Fully discrete problem Let us introduce the following spaces: V0h (Ωnh ) = {vh | vh ∈ Vh (Ωnh ), vh = 0 on Σh , vh |x1 =0 = vh |x1 =L } Z qh dx = 0, qh |x1 =0 = qh |x1 =L } P0h (Ωnh ) = {qh | qh ∈ Ph (Ωnh ), n+1/5,− Vh n+1/5,−
Γh
c n+3/5,+ W h b n+3/5,+ Γ h
Ωn h
= {vh | vh ∈
Vh (Ωnh ),
n+1/5,−
vh = 0 on Γh
n+1/5
= {x | x ∈ Γnh , uh
· n(x) < 0} n+3/5,+
b h ), vh = 0 on Γ b = {vh | vh ∈ Vh (Ω h n+3/5
bh , w = {x | x ∈ Γ h
, vh |x1 =0 = vh |x1 =L }
, vh |x1 =0 = vh |x1 =L },
· e2 > 0}.
Then, the fully discretized problem reads as follows. n+1/5
n+1/5
1. Find uh ∈ V0h (Ωnh ) and ph ∈ P0h (Ωnh ) such that Z Z Z n+1/5 uh − unh µ n+1/5 n+1/5 ̺ · v dx + ∇u : ∇v dx − ph ∇ · vh dx h h h n n n ∆t 2 Ωh Ωh Ωh Z L Z L n 2 −1/2 ∂1 ηhn (1 + |∂1 ηhn |2 )−1/2 ∂1 v2h dx1 (1 + |∂ η | ) ∂ v dx − s = −s 1 h 1 1h 1 0 0 Z Z µ n t + (∇uh ) : ∇vh dx, ∀vh ∈ V0h (Ωnh ) , g · vh dx − n n 2 Ω Ω h h Z n+1/5 qh ∇ · uh dx = 0, ∀qh ∈ Ph (Ωnh ) , Ωn h
where Remark 1 has been taken into account.
(5.1)
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
13
n+2/5
2. Find uh ∈ Vh (Ωnh ) via the solution of the following pure advection probn lem in Ωh × (tn , tn+1 ) Z Z n+1/5 ∂t uh · vh dx + (uh · ∇)uh · vh dx = 0, n n Ωh Ωh n+1/5 uh (tn ) = uh , n+1/5 n+1/5,− uh (t) = uh on Γh × (tn , tn+1 ) ,
n+1/5,−
∀vh ∈ Vh
, (5.2)
n+2/5
and then set uh = uh (tn+1 ). n+3/5 n+3/5 3. Find uh ∈ V0h (Ωnh ) and ph ∈ P0h (Ωnh ) such that Z Z Z n+3/5 uh − unh µ n+3/5 n+3/5 · vh dx + ∇uh : ∇vh dx − ph ∇ · vh dx ̺ n n n ∆t 2 Ω Ω Ω h h h Z µ n+2/5 t =− (∇uh ) : ∇vh dx, ∀vh ∈ V0h (Ωnh ) , n 2 Ωh Z n+3/5 qh ∇ · uh dx = 0, ∀q ∈ Ph (Ωnh ) . Ωn h
(5.3) b 4. Find ∈ Vh (Ωh ) via the solution of the following pure advection probb h × (tn , tn+1 ) lem in Ω n+4/5 ˆh u
Z Z n+3/5 b ˆh J(η)dξ = 0, ˆ ˆ (wh · ∇)ˆ uh · v u · v J(η)dξ − ∂ h t h b b Ωh Ωh n+3/5 ˆ h (tn ) = u ˆh u , n+3/5 b n+3/5,+ × (tn , tn+1 ) , ˆ h (t) = u ˆh u on Γ h
n+3/5,+
c ∀ˆ vh ∈ W h
,
(5.4) n+4/5 ˆh ˆ h (tn+1 ) and pn+4/5 = pn+1/5 + pn+1/3 . and then set u =u 5. Find the new position of the free surface ηhn+1 ∈ Sh by solving the following problem in (0, L) × (tn , tn+1 )
Z
0
L
∂t ηh φh dξ1 +
Z
n n ηh (ξ1 , t ) = ηh (ξ1 ),
0
L
n+4/5
(ˆ u1h
n+4/5
(ξ1 , H)∂1 ηh − u ˆ2h
(ξ1 , H)) φh dξ1 = 0,
∀φh ∈ Sh ,
(5.5)
14
Roland Glowinski and Giovanna Guidoboni
JMFM
b h trough the The new domain Ωn+1 is obtained from the reference domain Ω h mapping A−1 defined as: tn+1 ,h A−1 tn+1 ,h :
b h ∈ R2 Ω
→ Ωn+1 ∈ R2 h
x1 = ξ1 ξ = (ξ1 , ξ2 ) → x = (x1 , x2 ) = η n+1 (ξ1 ) x2 = h ξ2 . H
(5.6)
The velocity un+1 and the pressure pn+1 defined on the new domain Ωn+1 h h h are obtained as follows n+4/5
ˆh un+1 (x1 , x2 ) = u h
(x1 , ηhn (x1 )x2 /η n+1 (x1 )),
n+4/5
pn+1 (x1 , x2 ) = ph h
(x1 , ηhn (x1 )x2 /ηhn+1 ),
(5.7) (5.8)
where (x1 , x2 ) ∈ Ωn+1 . h The backward Euler’s method has been used to derive the time discretization in (5.1) and (5.3). These problem can be seen as a ”generalized” discrete Stokes problem for which efficient solution methods already exist as shown in, e.g., [8]. The advection step (3.2) has been solved using a wave-like equation method [8, 10, 17]: this approach preserves the hyperbolic nature of advection, it does not introduce numerical dissipation and it is easy to implement. It consists in replacing problem (5.2) by the following wave-like equation one in Ωnh × (tn , tn+1 ): Z Z n+1/5 n+1/5 2 ∂ u · v dx + (uh · ∇)uh · (uh · ∇)vh dx h t h n n Ωh Ωh Z n+1/5 n+1/5,− + (uh · n)(∂t uh · vh )dΓnh = 0, ∀vh ∈ V0h , n+1/5,− Γn \Γ h h n+1/5 uh (tn ) = uh , Z Z n+1/5 n+1/5,− n ∂ u (t ) · v dx = − (un+1/5 · ∇)uh · vh dx, ∀vh ∈ V0h , t h h n n Ωh Ωh n+1/5,− n+1/5 n+1/5,− ∂t uh (tn ) ∈ V0h , uh = uh on Γh × (tn , tn+1 ). (5.9) We have solved problem (5.9) using a second order accurate time discretization scheme which is discussed in, e.g., [8], Chapter 6, and [17]. To solve the trasport problem (5.4) we use again the wave equation approach. ˆ and therefore the We observe that (5.4) does not contain x1 differentiation of u problem reduces to the solution of a family (infinite for the continuous problem, finite for the discrete ones) of transport problems in one space dimension along ˆ is solution of a the vertical direction. Then for ξ1 ∈ [0, L), each component of u
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
transport problem of the following form: ∂ϕ ∂ϕ = 0 on (0, L) × (tn , tn+1 ), − aξ2 ∂t ∂ξ2 ϕ(tn ) = ϕ0 , ϕ(H, t) = b if a > 0, t ∈ (tn , tn+1 ),
15
(5.10)
with a and b constant with respect to ξ2 and t. We observe that any smooth solution of (5.10) verifies: 2 ∂ ϕ ∂ϕ ∂ n n+1 2 ), ∂t2 − a ξ2 ξ2 ξ2 ξ2 = 0 on (0, L) × (t , t ϕ(tn ) = ϕ0 , ∂ϕ ∂ϕ0 , (tn ) = aξ2 (5.11) ∂t ∂ξ2 n n+1 ϕ(H, t) = b if a > 0, t ∈ (t , t ), ∂ϕ ∂ϕ = 0 if a ≤ 0, t ∈ (tn , tn+1 ). aξ2 ∂t − aξ2 ∂ξ 2 ξ2 =H
From a practical point of view, we take advantage of the following variational formulation of (5.10), better suited for finite element treatment. 1. Case a > 0. Z H 2 Z H ∂ ϕ v ∂ϕ ∂v 2 dξ + a dξ2 = 0, ξ2 2 2 ξ ∂t ∂ξ 2 2 ∂ξ2 0 0 ϕ(t0 ) = ϕ0 ,
∀v ∈ D(0, H), a.e. on (t0 , tf ),
∂ϕ0 ∂ϕ , (t0 ) = aξ2 ∂t ∂ξ2 ϕ(H, t) = f,
with D(0, H) = {v | v ∈ C ∞ [0, H], v has compact support in (0, H) }. 2. Case a < 0. Z H 2 Z H ∂ϕ ∂ϕ ∂v ∂ ϕ v 2 dξ2 + a dξ2 − a (H, t)v(H) = 0, ξ2 2 ξ ∂t ∂ξ ∂ξ ∂t 2 2 2 0 0 ∀v ∈ D(0, H), a.e. on (t0 , tf ), ϕ(t0 ) = ϕ0 , ∂ϕ0 ∂ϕ , (t0 ) = aξ2 ∂t ∂ξ2
with V0 = {v | v ∈ C ∞ [0, H], v has compact support in (0, H] }.
(5.12)
(5.13)
16
Roland Glowinski and Giovanna Guidoboni
JMFM
Finally, problem (5.5) is solved using a second order Taylor-Galerkin discretization on (0, L) × (tn , tn+1 ), namely Z L Z L n+4/5 n+4/5 (ηhn+1 − ηhn ) φh dx1 = ∆t (u2h − u1h ∂1 ηhn ) φh dx1 0 0 Z (∆t)2 L n+4/5 n+4/5 (5.14) u1h ∂1 u2h φh dx1 − 2 0 Z (∆t)2 L n+4/5 n+4/5 n+4/5 − u1h ∂1 ηhn (∂1 u1h φh + u1h ∂1 φh ) dx1 , 2 0 where the velocity is computed at Γnh . Remark 2. For a general approach concerning the approximation of free surface problems for incompressible viscous fluids with surface tension see, e.g., [4].
6. Numerical experiments Problem (2.4)-(2.10) always admits a steady solution. Indeed, assuming that u = U (x2 )e1 , with e1 = (1, 0), p = P (x2 ) and η(x1 , t) ≡ H, we get U=
̺g sin β (2Hx2 − x22 ), 2µ
P = −̺g cos β(x2 − H) in (0, L) × (0, H). (6.1)
Nishida et al. in [15] proved the existence of time periodic solutions bifurcating from the basic flow (6.1) for Reynolds numbers large enough. The Reynolds number Re is defined as Re = HUmax ̺/µ, where H is the average height of the fluid layer and Umax is the velocity of the basic flow at the free surface, obtained from (6.1) with x2 = H, namely H 2 ̺g sin β Umax = , (6.2) 2µ implying that H 3 ̺2 g sin β Re = . (6.3) 2µ2 Our goal is to explore the features of the bifurcating solutions using the numerical algorithm described in the previous sections. For the time discretization we have taken ∆t = 2.5 · 10−4 . As mentioned in Section 5, we solve the advection steps (5.2) and (5.4) using a wave-like equation method which requires the specification of time substeps: here we have chosen both of them equal to ∆t/5. In our numerical experiments, the dimensions of the reference rectangular domain are H = 0.15 and L = 3, the density of the fluid is ̺ = 1, and the angle of inclination is β = 45◦ . By direct simulation, we estimated the critical Reynolds number Rec as a function of the surface tension. In particular, we considered s = 1, 2, 3. The results are reported in Table 1. It is clear that Rec is an increasing function of s, which means that the surface tension has a stabilizing effect, as shown in [13].
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
17
range of Rec range for µc surface tension s s = 1 1.68 < Rec < 1.98 1/12 < µc < 1/13 s = 2 2.29 < Rec < 2.63 1/14 < µc < 1/15 3.00 < Rec < 3.38 1/16 < µc < 1/17 s = 3 Table 1. Range of the critical Reynolds number Rec and the correspondent viscosities for several values of the surface tension s.
The system reaches a steady state whenever the Reynolds number is below the critical one. In the simulations, the solution is considered a steady state when the following condition is satisfied: kun − un−1 kL2 (Ωn ) < 10−7 . kun kL2 (Ωn )
(6.4)
In Figure 5 we compare the analytical velocity profile given by (6.1) and the one obtained numerically with our algorithm (here s = 1, µ = 1/11, and Re = 1.42). The agreement between the two profiles is indeed very satisfactory. When the Reynolds number exceeds the critical threshold, periodic solutions establish. In Figure 8 this phenomenon is tracked by the vertical oscillations of the middle point on the free surface (here s = 1, µ = 1/15, and Re = 2.63). In Figure 9 we have plotted five pictures of the streamlines contours showing the evolution of the solution at Re = 2.63 and s = 1 during one whole period. The wavelength of the waves on the free surface is very long at the early stages of the onset of the instability, as shown in Figure 6. As time goes by, the system enters a periodic regime and the surface waves develop into a profile whose shape, amplitude and frequency depend on the Reynolds number and on the surface tension. Numerical results for s = 1 and different Reynolds numbers are reported in Table 2, and they show that both the amplitude and the frequency of the waves are increasing functions of Re. Moreover, the shape of the waves steepen as the Reynolds number increases, as shown in Figure 7. We repeated the computations for s = 1 and Re = 1.98 using different mesh sizes and time steps. The results are reported in Table 3 and they show a very good agreement in terms of amplitude, frequency and period of the periodic oscillations. Animations of the numerical simulations described in this article are available at http://www.math.uh.edu/∼gio/
References [1] M. Behr, F. Abraham, Free-Surface Flow Simulations in the Presence of Inclined Walls, Comput. Methods Appl. Mech. and Engrg. 191 (2002) 5467-5483 [2] T.B. Benjamin Wave formation in laminar flow down an inclined plane, J. Fluid Mech. 2 (1957) 554-574
18
Roland Glowinski and Giovanna Guidoboni
JMFM
0.14
0.12
x2
0.1
0.08
0.06
0.04
0.02
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
Figure 5. Comparison between the velocity profile at the steady state computed analytically (solid line) and numerically (×). Here s = 1, µ = 1/11, and Re = 1.42. Re amplitude period frequency 1.98 0.03215 1.63150 3.85117 2.63 0.06956 1.41775 4.43180 0.13927 1.15650 5.43293 3.79 Table 2. Amplitude and frequency for the oscillations of the middle point on the free surface for s = 1 and several Reynolds number beyond the Hopf bifurcation. amplitude period frequency (hp , ∆t) (1/100, 2.5 × 10−4 ) 0.03215 1.63150 3.85117 (1/64, 2.5 × 10−4 ) 0.03210 1.63275 3.84822 0.03176 1.6291 3.85691 (1/100, 1.25 × 10−4 ) Table 3. Amplitude and frequency for the oscillations of the middle point on the free surface for s = 1, µ = 1/13 and Re = 1.98 for different combinations of mesh sizes and time steps.
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
19
0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0
0
0.5
1
1.5
2
2.5
3
Figure 6. The instability begins with waves characterized by very long wavelengths. Here s = 1, µ = 1/15, and Re = 2.63. [3] M. Bercovier, O. Pironneau, Error estimates for finite element method solution of the Stokes problem in primitive variables, Numer. Math. 33 (1979) 211-224 [4] A. Caboussat, A numerical method for the simulation of free surface flows with surface tension, Computers and Fluids 35-10 (2006) 1205-1216 [5] R.W. Chin, F.H. Abernathy, R. Bertschy, Gravity and shear wave stability of free surface flows. Part1. Numerical calculations, J. Fluid Mech. 168 (1986) 501-513 [6] W. Dettmer, D. Peri´c, A computational framework for free surface accounting for surface tension, Comput. Methods Appl. Mech. Engrg. 195 (2006) 3038-3071 [7] J. Donea, Arbitrary Lagrangian Eulerian Methods, in Computational Methods for Transient Analysis, Vol. I of Computational Methods in Mechanics. North-Holland, Elsevier (1983) [8] R. Glowinski, Finite Element Methods for Incompressible Viscous Flow, in Handbook of Numerical Analysis, Vol. IX, P.G. Ciarlet and J.-L. Lions eds., North-Holland, Amsterdam (2003) [9] R. Glowinski, L.H. Ju´ arez, Finite element method and operator splitting for a timedependent viscous incompressible free surface flow, Comput. Fluid Dyn. J. 12-3 (2003) 459-468 [10] R. Glowinski, G. Guidoboni, T.W. Pan, Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity, J. Comput. Phys. 216-1 (2006) 76-91
20
Roland Glowinski and Giovanna Guidoboni
JMFM
0.24
0.22
0.2
0.18
0.16
0.14
0.12
0.1
0.08
0
0.5
1
1.5
2
2.5
3
Figure 7. Shape of the waves on the free surface in the periodic regime for Re = 1.98 and µ = 1/13 (solid line), Re = 2.63 and µ = 1/15 (dashed line), Re = 3.79 and µ = 1/18 (dashed-dotted line), when s = 1. [11] D.Rh. Gwynllyw, D.H. Peregrine, Numerical Simulation of Stokes flow down an inclined plane, Proc. R. Soc. Lond. A 452-1946 (1996) 543-565 [12] T.J. Hughes, W.K. Liu, T.K. Zimmermann, Lagrangian-eulerian finite element formulation for incompressible viscous flows, Comp. Meth. Appl. Mech. Engrg. 29 (1981) 329-349 [13] J. Liu, J.D. Paul, J.P. Gollub, Measurements of the primary instabilites of film flows, J. Fluid Mech. 250 (1993) 69-101 [14] G.I. Marchuk, Splitting and alternate direction methods, in Handbook of Numerical Analysis, Vol. I, P.G. Ciarlet and J.-L. Lions eds., North-Holland, Amsterdam (1990) [15] T. Nishida, Y. Teramoto, H. Yoshihara, Hopf bifurcation in viscous incompressible flow down an inclined plane, J. Math. Fluid Mech. 7 (2005) 29-71 [16] M. Padula, V.A. Solonnikov, On Rayleigh-Taylor stability, Ann. Univ. Ferrara Sez. VII (N.S.) 46 (2000), 307-336 [17] T.W. Pan and R. Glowinski, A projection/wave-like equation method for the numerical simulation of incompressible viscous fluid flow modeled by the Navier-Stokes equations , Comput. Fluid Dyn. J. 9-2 (2000) 28-42 [18] O. Pironneau, Finite Element Methods for Fluids, J. Wiley, Chichester (1989)
Vol. 99 (9999)
Hopf bifurcation for a free surface flow
21
0.18
0.17
0.16
0.15
0.14
0.13
0.12 0
5
10
15
20
25
30
35
Figure 8. Variation in time of height of the middle point of the free surface (s = 1, µ = 1/15 and Re = 2.63).
Figure 9. Streamlines of the solution during one period, obtained for s = 1, µ = 1/15 and Re = 2.63. Acknowledgment The authors would like to thank Professors H. Juarez, M. Padula, T.W. Pan, M. Picasso and J. Rappaz for helpful comments and suggestions. The support of
22
Roland Glowinski and Giovanna Guidoboni
JMFM
NSF (Grant DMS-9973318, 0209066, and CCR-9902035) and DOE/LACSI (Grant R71700K-292-000-99) is also acknowledged. Roland Glowinski Department of Mathematics, University of Houston, Houston, Texas 77204-3476, USA and Laboratoire Jacques-Louis Lions, Universit´e P. et M. Curie, 4 Place Jussieu, 75005, Paris, France e-mail: [email protected] Giovanna Guidoboni Department of Mathematics, University of Houston, Houston, Texas 77204-3476, USA e-mail: [email protected]
