- Email: [email protected]
Simulation of an oil slick movement using a shallow water model
Mathematics and Computers in Simulation 76 (2007) 155–160
Simulation of an oil slick movement using a shallow water model B. Di Martino, M. Peybernes ∗ UMR CNRS 6134, Universit´e de Corse, 20250 Corte, France Available online 23 January 2007
Abstract The behaviour of an oil slick on the sea surface is studied using a shallow water model with free boundary. This model is solved by coupling the Arbitrary Lagrangian Eulerian method (ALE) with the characteristic method and by using a Galerkin spatial discretization. A special basis, which permits to obtain easily the solution by solving a scalar eigenvalue problem, is then proposed. Finally, a numerical experiment is presented in a realistic configuration. © 2007 IMAC. Published by Elsevier B.V. All rights reserved. AMS Classification: 35Q30; 65M60; 76D03 Keywords: Shallow water equations; Free boundary problem; Galerkin and ALE discretization
1. Introduction This paper presents the modelization and the simulation of an oil slick movement on the sea surface. The objective is to predict the movement of pollutants generated by marine spills and to have a good representation of the surface occupied by these pollutants. The oil slick is modelized with a shallow water model written under the following form ∂u + (u · ∇)u + g∇h − νu = f1 , ∂t
(1)
∂h + div(uh) = f2 . ∂t
(2)
where u is the mean horizontal velocity (integrated on the vertical) of the oil slick, h the thickness of the oil slick, g the gravity coefficient and ν is the viscosity coefficient. f1 and f2 are, respectively, the external forcing (for example, wind stress, sea surface current stress, sea surface inclination) and a source or a vanishing of oil (for example, by evaporation or dissolution). This system is studied in a two-dimensional time dependent domain t , with a boundary t . Physically, this domain represents the surface occupied by the oil slick at time t. The displacement d of the boundary t is given by the relation ∂d = u on t . ∂t ∗
Corresponding author. Tel.: +33 4 95 45 00 41. E-mail address: [email protected] (M. Peybernes).
0378-4754/$32.00 © 2007 IMAC. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2007.01.022
156
B. Di Martino, M. Peybernes / Mathematics and Computers in Simulation 76 (2007) 155–160
We complete this boundary condition by a condition on the divergence and the curl of the velocity on t . In this approach, we have chosen to use h = 0;
curl u = 0;
div u = 0 on t .
The paper is organized as follows. In Section 2, it is show how the Arbitrary Lagrangian Eulerian (ALE) method can be applied for the previous system. In Section 3, a numerical experiment with a realistic condition is reported. Comments and conclusions are contained in Section 4. 2. The numerical method 2.1. Review of the ALE method In order to avoid unnecessary notations, we do not take into account the terms f1 and f2 in this section. The ALE method has been mainly used for the resolution of the Navier-Stokes equations in a moving domain [4]. We assume that at the initial time t0 , the fluid domain is overlapped by a regular mesh (for example, an usual finite element mesh). We assume that on the boundary, each point of this mesh has the velocity of the fluid. The velocity of the internal points can be chosen arbitrarily, the unique condition is the regularity of the mesh at each time. Then, for any time t ∈ [0, T ], the velocity ct of the moving domain can be defined by solving the following problem: ct = 0 in t
(3)
ct = u(x, t) on t
(4)
and we consider the following mapping ˜ → R2 c:Q (x, t) → c(x, t) = ct (x), ˜ = ∪t ∈ [0,T ] (t × {t}). The relation between the different domains ti is given by the mapping where Q C(·, t1 , t2 ) : t1 → t2 ,
x1 → x2 = C(x1 , t1 , t2 )
where C(·, t1 , t2 ) is the characteristic curve from (x1 , t1 ) to (x2 , t2 ) corresponding to the velocity ct in the space-time ˜ For each time τ, we denote by uτ and hτ the ALE velocity and the ALE thickness: domain Q. uτ (x, t) = u(C(x, τ, t), t),
hτ (x, t) = h(C(x, τ, t), t).
(5)
With these notations, we can write the ALE formulation of the problem ∂uτ + (uτ − cτ )∇uτ − μuτ + g∇hτ = 0, ∂t
(6)
∂hτ + (uτ − cτ )∇hτ + hτ div uτ = 0. ∂t
(7)
We note t the time step, n = t n with boundary n and n (x) = g m (x, t n ) for x ∈ m gm t gn (x) = gnn (x) = g(x, t n ) for x ∈ n
In order to ensure the positivity of ht n , we renormalize the continuity equation as follows ∂ log ht n + (ut n − ct n )∇ log ht n + div ut n = 0. ∂t
(8)
We denote by U n (respectively, H n and Cn ) the approximation of the exact solution ut n (respectively, ht n and ct n ).
B. Di Martino, M. Peybernes / Mathematics and Computers in Simulation 76 (2007) 155–160
157
˜ n (x) = U n (Xn (x, t n )) and H ˜ n (x) = H n (Xn (x, t n )) with Xn (x, ·) the characteristic curve solution Then, we note U of:
⎧ n ⎨ ∂X (x, t) = (U n − Cn )(Xn (x, t)) ∂t ⎩ n X (x, t n+1 ) = x.
In our study, we approximate the foot of the characteristic by Xn (x, t n ) x − (U n − Cn )(x)t. Approximating the Lagrangian derivative by a first-order Euler scheme, we obtain: ˜ n (x) Unn+1 (x) − νtUnn+1 (x) + gt ∇Hnn+1 (x) = U ˜ n (x) log Hnn+1 (x) + t div Unn+1 (x) = log H div Unn+1 = curl Unn+1 = Hnn+1 = 0
in n
in n
in n
(9) (10) (11)
2.2. Spatial discretization The spatial discretization of the previous problem is based on the Galerkin method [5]. For that, we use the following orthogonal decomposition L2 (n )2 = ∇H01 (n ) ⊕ curl H01 (n ) ⊕ ∇H(n ) where H is the intersection of H 1 and the space of harmonic functions. Thus, Unn+1 can be decomposed as follows + curl qnn+1 + ∇rnn+1 Unn+1 = ∇pn+1 n n+1 ∈ H 1 (n ) and r n+1 ∈ H. Since, div U n+1 = curl U n+1 = H n+1 = 0 on n , we deduce that pn+1 with pn+1 n , qn n n n n n 0 n+1 and qn vanish on the boundary. The problem is then: ⎧ ˜ n in n ∇(pn+1 − ν t pn+1 + g t ∇Hnn+1 ) = P∇H 1 U ⎪ n n ⎪ 0 ⎨ ˜n curl (qnn+1 − νtqnn+1 ) = Pcurl H 1 U in n 0 ⎪ ⎪ ⎩ ∇r n+1 = P U ˜ n − P 1U ˜n−P ˜n in n 1U ∇H
n
∇H0
curl H0
with the boundary conditions Hnn+1 = pn+1 = qnn+1 = 0 on n . With this decomposition, we have directly the projecn n+1 n ˜ . The two first equations can be replaced by tion of Un on ∇H according to U − νtpn+1 + gtHnn+1 = θ in n pn+1 n n qnn+1 − νt qnn+1 = δ
in n
˜ n on ∇H 1 and curl H 1 . where θ and δ come from the projection of U 0 0 The previous system is solved by using the Galerkin method with pn+1 n
=
Np
αni ei ,
qnn+1
=
Nq
i=1
βin ei
i=1
where {ek }k ∈ N is a basis of verifying the eigenvalue problem: −ek = λk ek in n H01
ek = 0
on n
Then, we approximate Unn+1 by Unn+1 =
Np i=1
αni ∇ei +
Nq i=1
βin curl ei + ∇rnn+1
(12)
158
B. Di Martino, M. Peybernes / Mathematics and Computers in Simulation 76 (2007) 155–160
and Hnn+1 by a characteristic method: ˜ nn+1 e−tdiv Unn+1 Hnn+1 = H
(13)
3. Numerical experiment The previous method is tested in order to simulate the movement of an oil slick in the dam of Calacuccia (Corsica). We assume that the sea surface velocity is known (numerous methods to compute this can be found in refs. [1,2]). The oil slick moves under three different forcings: the wind stress under the form Cw (u − W)|u − W| (where W is the wind velocity), the sea stress under the form Cs (u − us )|u − us | (where us is the sea velocity on the surface of the sea) and the gravity effects. At initial time t0 , oil slick is assumed to be circular (Fig. 1). After 3 h, the boundary of the oil has moved under the different forcings (Fig. 2). The thickness of the oil slick is smaller and the surface occupied by the oil is bigger. The simulation is stopped after 5 h when the oil reaches the boundary of the lake Calacuccia (Fig. 3). The resulting mesh after 3 h is presented in the Fig. 4.
Fig. 1. Initial location of the pollutant.
Fig. 2. Position after 3 h.
B. Di Martino, M. Peybernes / Mathematics and Computers in Simulation 76 (2007) 155–160
Fig. 3. Position after 5 h.
Fig. 4. Slick mesh after 3 h.
159
160
B. Di Martino, M. Peybernes / Mathematics and Computers in Simulation 76 (2007) 155–160
4. Conclusions In this study, we have applied a coupled Galerkin–ALE method to a shallow water problem on a domain with free boundary. The choice of these two methods permits to obtain a good representation of the oil slick movement on the surface of the sea. But the approach presented here uses a very simple boundary condition, especially by taking div u = 0 and h = 0 on t . More generally, the good condition is a condition on the trace of the stress tensor on the boundary given by γ(h − νdivu) = A(d), where γ is the trace operator and A is an appropriate operator taking into account the surface stress between sea and oil on the boundary t . This difficulty is going to be studied in a further work. Another difficulty can be observed when the deformation of the oil slick gives a very complex domain and we have not a sufficiently regular mesh (this is the case when a part of the oil reaches to the coast). Computation needs then to be stopped and we have not find at this time a good remedy. Acknowledgments This work was supported in part by the European Community (Interreg IIIA scheme) and an earmarked studentship from collectivit´e Territoriale de Corse. References [1] F. Bosseur, B. Di Martino, P. Orenga, Resolution by Galerkin method with a special basis of a geophysical flow in open sea: a Calvi’s bay simulation, Appl. Math. Model. 24 (2000) 73–94. [2] F. Flori, C. Giacomoni, P. Orenga, Convergence of a Characteristic-Galerkin scheme for a shallow water problem, Math. Comput. Model., in press. [4] B. Maury, Characteristics ALE method for the Unsteady 3D Navier Stokes Equations with a free surface, Int. J. Comput. Fluid Dyn. 6 (1996) 175–188. [5] P. Orenga, Construction d’une base sp´eciale pour la r´esolution de quelques probl`emes d’oc´eanographie physique en dimension deux, CRAS, vol. 314, Springer-Verlag, 1996, pp. 587–590.
Simulation of an oil slick movement using a shallow water model B. Di Martino, M. Peybernes ∗ UMR CNRS 6134, Universit´e de Corse, 20250 Corte, France Available online 23 January 2007
Abstract The behaviour of an oil slick on the sea surface is studied using a shallow water model with free boundary. This model is solved by coupling the Arbitrary Lagrangian Eulerian method (ALE) with the characteristic method and by using a Galerkin spatial discretization. A special basis, which permits to obtain easily the solution by solving a scalar eigenvalue problem, is then proposed. Finally, a numerical experiment is presented in a realistic configuration. © 2007 IMAC. Published by Elsevier B.V. All rights reserved. AMS Classification: 35Q30; 65M60; 76D03 Keywords: Shallow water equations; Free boundary problem; Galerkin and ALE discretization
1. Introduction This paper presents the modelization and the simulation of an oil slick movement on the sea surface. The objective is to predict the movement of pollutants generated by marine spills and to have a good representation of the surface occupied by these pollutants. The oil slick is modelized with a shallow water model written under the following form ∂u + (u · ∇)u + g∇h − νu = f1 , ∂t
(1)
∂h + div(uh) = f2 . ∂t
(2)
where u is the mean horizontal velocity (integrated on the vertical) of the oil slick, h the thickness of the oil slick, g the gravity coefficient and ν is the viscosity coefficient. f1 and f2 are, respectively, the external forcing (for example, wind stress, sea surface current stress, sea surface inclination) and a source or a vanishing of oil (for example, by evaporation or dissolution). This system is studied in a two-dimensional time dependent domain t , with a boundary t . Physically, this domain represents the surface occupied by the oil slick at time t. The displacement d of the boundary t is given by the relation ∂d = u on t . ∂t ∗
Corresponding author. Tel.: +33 4 95 45 00 41. E-mail address: [email protected] (M. Peybernes).
0378-4754/$32.00 © 2007 IMAC. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2007.01.022
156
B. Di Martino, M. Peybernes / Mathematics and Computers in Simulation 76 (2007) 155–160
We complete this boundary condition by a condition on the divergence and the curl of the velocity on t . In this approach, we have chosen to use h = 0;
curl u = 0;
div u = 0 on t .
The paper is organized as follows. In Section 2, it is show how the Arbitrary Lagrangian Eulerian (ALE) method can be applied for the previous system. In Section 3, a numerical experiment with a realistic condition is reported. Comments and conclusions are contained in Section 4. 2. The numerical method 2.1. Review of the ALE method In order to avoid unnecessary notations, we do not take into account the terms f1 and f2 in this section. The ALE method has been mainly used for the resolution of the Navier-Stokes equations in a moving domain [4]. We assume that at the initial time t0 , the fluid domain is overlapped by a regular mesh (for example, an usual finite element mesh). We assume that on the boundary, each point of this mesh has the velocity of the fluid. The velocity of the internal points can be chosen arbitrarily, the unique condition is the regularity of the mesh at each time. Then, for any time t ∈ [0, T ], the velocity ct of the moving domain can be defined by solving the following problem: ct = 0 in t
(3)
ct = u(x, t) on t
(4)
and we consider the following mapping ˜ → R2 c:Q (x, t) → c(x, t) = ct (x), ˜ = ∪t ∈ [0,T ] (t × {t}). The relation between the different domains ti is given by the mapping where Q C(·, t1 , t2 ) : t1 → t2 ,
x1 → x2 = C(x1 , t1 , t2 )
where C(·, t1 , t2 ) is the characteristic curve from (x1 , t1 ) to (x2 , t2 ) corresponding to the velocity ct in the space-time ˜ For each time τ, we denote by uτ and hτ the ALE velocity and the ALE thickness: domain Q. uτ (x, t) = u(C(x, τ, t), t),
hτ (x, t) = h(C(x, τ, t), t).
(5)
With these notations, we can write the ALE formulation of the problem ∂uτ + (uτ − cτ )∇uτ − μuτ + g∇hτ = 0, ∂t
(6)
∂hτ + (uτ − cτ )∇hτ + hτ div uτ = 0. ∂t
(7)
We note t the time step, n = t n with boundary n and n (x) = g m (x, t n ) for x ∈ m gm t gn (x) = gnn (x) = g(x, t n ) for x ∈ n
In order to ensure the positivity of ht n , we renormalize the continuity equation as follows ∂ log ht n + (ut n − ct n )∇ log ht n + div ut n = 0. ∂t
(8)
We denote by U n (respectively, H n and Cn ) the approximation of the exact solution ut n (respectively, ht n and ct n ).
B. Di Martino, M. Peybernes / Mathematics and Computers in Simulation 76 (2007) 155–160
157
˜ n (x) = U n (Xn (x, t n )) and H ˜ n (x) = H n (Xn (x, t n )) with Xn (x, ·) the characteristic curve solution Then, we note U of:
⎧ n ⎨ ∂X (x, t) = (U n − Cn )(Xn (x, t)) ∂t ⎩ n X (x, t n+1 ) = x.
In our study, we approximate the foot of the characteristic by Xn (x, t n ) x − (U n − Cn )(x)t. Approximating the Lagrangian derivative by a first-order Euler scheme, we obtain: ˜ n (x) Unn+1 (x) − νtUnn+1 (x) + gt ∇Hnn+1 (x) = U ˜ n (x) log Hnn+1 (x) + t div Unn+1 (x) = log H div Unn+1 = curl Unn+1 = Hnn+1 = 0
in n
in n
in n
(9) (10) (11)
2.2. Spatial discretization The spatial discretization of the previous problem is based on the Galerkin method [5]. For that, we use the following orthogonal decomposition L2 (n )2 = ∇H01 (n ) ⊕ curl H01 (n ) ⊕ ∇H(n ) where H is the intersection of H 1 and the space of harmonic functions. Thus, Unn+1 can be decomposed as follows + curl qnn+1 + ∇rnn+1 Unn+1 = ∇pn+1 n n+1 ∈ H 1 (n ) and r n+1 ∈ H. Since, div U n+1 = curl U n+1 = H n+1 = 0 on n , we deduce that pn+1 with pn+1 n , qn n n n n n 0 n+1 and qn vanish on the boundary. The problem is then: ⎧ ˜ n in n ∇(pn+1 − ν t pn+1 + g t ∇Hnn+1 ) = P∇H 1 U ⎪ n n ⎪ 0 ⎨ ˜n curl (qnn+1 − νtqnn+1 ) = Pcurl H 1 U in n 0 ⎪ ⎪ ⎩ ∇r n+1 = P U ˜ n − P 1U ˜n−P ˜n in n 1U ∇H
n
∇H0
curl H0
with the boundary conditions Hnn+1 = pn+1 = qnn+1 = 0 on n . With this decomposition, we have directly the projecn n+1 n ˜ . The two first equations can be replaced by tion of Un on ∇H according to U − νtpn+1 + gtHnn+1 = θ in n pn+1 n n qnn+1 − νt qnn+1 = δ
in n
˜ n on ∇H 1 and curl H 1 . where θ and δ come from the projection of U 0 0 The previous system is solved by using the Galerkin method with pn+1 n
=
Np
αni ei ,
qnn+1
=
Nq
i=1
βin ei
i=1
where {ek }k ∈ N is a basis of verifying the eigenvalue problem: −ek = λk ek in n H01
ek = 0
on n
Then, we approximate Unn+1 by Unn+1 =
Np i=1
αni ∇ei +
Nq i=1
βin curl ei + ∇rnn+1
(12)
158
B. Di Martino, M. Peybernes / Mathematics and Computers in Simulation 76 (2007) 155–160
and Hnn+1 by a characteristic method: ˜ nn+1 e−tdiv Unn+1 Hnn+1 = H
(13)
3. Numerical experiment The previous method is tested in order to simulate the movement of an oil slick in the dam of Calacuccia (Corsica). We assume that the sea surface velocity is known (numerous methods to compute this can be found in refs. [1,2]). The oil slick moves under three different forcings: the wind stress under the form Cw (u − W)|u − W| (where W is the wind velocity), the sea stress under the form Cs (u − us )|u − us | (where us is the sea velocity on the surface of the sea) and the gravity effects. At initial time t0 , oil slick is assumed to be circular (Fig. 1). After 3 h, the boundary of the oil has moved under the different forcings (Fig. 2). The thickness of the oil slick is smaller and the surface occupied by the oil is bigger. The simulation is stopped after 5 h when the oil reaches the boundary of the lake Calacuccia (Fig. 3). The resulting mesh after 3 h is presented in the Fig. 4.
Fig. 1. Initial location of the pollutant.
Fig. 2. Position after 3 h.
B. Di Martino, M. Peybernes / Mathematics and Computers in Simulation 76 (2007) 155–160
Fig. 3. Position after 5 h.
Fig. 4. Slick mesh after 3 h.
159
160
B. Di Martino, M. Peybernes / Mathematics and Computers in Simulation 76 (2007) 155–160
4. Conclusions In this study, we have applied a coupled Galerkin–ALE method to a shallow water problem on a domain with free boundary. The choice of these two methods permits to obtain a good representation of the oil slick movement on the surface of the sea. But the approach presented here uses a very simple boundary condition, especially by taking div u = 0 and h = 0 on t . More generally, the good condition is a condition on the trace of the stress tensor on the boundary given by γ(h − νdivu) = A(d), where γ is the trace operator and A is an appropriate operator taking into account the surface stress between sea and oil on the boundary t . This difficulty is going to be studied in a further work. Another difficulty can be observed when the deformation of the oil slick gives a very complex domain and we have not a sufficiently regular mesh (this is the case when a part of the oil reaches to the coast). Computation needs then to be stopped and we have not find at this time a good remedy. Acknowledgments This work was supported in part by the European Community (Interreg IIIA scheme) and an earmarked studentship from collectivit´e Territoriale de Corse. References [1] F. Bosseur, B. Di Martino, P. Orenga, Resolution by Galerkin method with a special basis of a geophysical flow in open sea: a Calvi’s bay simulation, Appl. Math. Model. 24 (2000) 73–94. [2] F. Flori, C. Giacomoni, P. Orenga, Convergence of a Characteristic-Galerkin scheme for a shallow water problem, Math. Comput. Model., in press. [4] B. Maury, Characteristics ALE method for the Unsteady 3D Navier Stokes Equations with a free surface, Int. J. Comput. Fluid Dyn. 6 (1996) 175–188. [5] P. Orenga, Construction d’une base sp´eciale pour la r´esolution de quelques probl`emes d’oc´eanographie physique en dimension deux, CRAS, vol. 314, Springer-Verlag, 1996, pp. 587–590.