CC BY
6
Вычислительные технологии
Том 24, № 1, 2019
Asymptotic behaviour of dam break flow for small times
D. I. Demirel1, A. Iafrati2, A. Korobkin3’!, O. Yilmaz1 1Izmir Institute of Technology, Dept. of Mathematics, 35430 Urla, Izmir, Turkey 2National Research Council of Italy Institute of Marine Engineering, 00128 Roma, Italy, 3University of East Anglia, School of Mathematics, NR4 7TJ Norwich, UK tCorresponding author: Korobkin, Alexander, e-mail: A.Korobkin@uea.ac.uk
Two dimensional impulsive flow of a fluid is studied within the potential flow theory. Initially the fluid is at rest and is held on one side of a vertical plate. The plate is withdrawn suddenly and gravity driven flow of the fluid starts. Attention is paid to the singular behaviour of the velocity field at the bottom point, where the vertical free surface meets the rigid bottom. The linear problem is solved by the Fourier series method. An inner region solution is found using Mellin transform at the bottom point.
The jet formation is observed at the bottom point. Also the discontinuity at the upper corner point is dealt with Lagrangian variables. For the second order outer problem, domain decomposition method is used. Comparison of the shapes of the free surfaces near the upper corner point with leading and second order solutions shows that the second order outer solution outer makes a larger difference in the vertical free surface than in the horizontal portion, compared with leading order solution.The complete picture of the shapes of the free surfaces using Lagrangian description for the upper part and Eulerian description for the bottom part at the second order is obtained.
Keywords: dam-break flow, free-surface flow, matched asymptotic expansions, domain decomposition method.
Cite: Demirel, D.I., Iafrati, A., Korobkin, A., Yilmaz, O. Asymptotic behaviour of dam break flow for small times // Computational Technologies. 2019. Vol. 24, No. 1.
P. 7-27. DOI: 10.25743/ICT.2019.24.1.002.
Introduction
The initial stages of the dam break flow, which is caused when a vertical dam at x' = 0, — H<y' < 0, holding the liquid at the semi infinite strip, x' > 0, —H < у' < 0, suddenly disappears, is investigated. At t' = 0, the fluid is at rest above a rigid bed with depth H. A cartesian coordinate system (х', у') with an origin at the free surface and positive x'-axis directed along the free surface is chosen as shown in Fig. 1. The resulting flow is gravity driven, two dimensional and potential. The liquid is assumed to be inviscid and incompressible.
There are two free-surfaces of the flow region, which vary in time and have to be determined as a part of solution. The upper horizontal part of the free-surface is denoted as у' = у'(x',t'), x' > 0 and the vertical part as x' = £'(y',t'). The flow region is bounded by these free-surfaces and by the rigid bottom у' = —H.
© ICT SB RAS, 2019
Free surface
_^ v'
Dam H
Fig. 1. Flow region at the initial time instant t' = 0
In this dam-break problem, we aim to construct a uniformly valid small-time solution by using matched asymptotic expansions. The solution in time as power series should be considered 'outer' solution, which should be corrected with 'inner' solution near the intersection point, where there is a logarithmic singularity of the free surface shape and the horizontal fluid velocity. The outer and inner asymptotic solutions have to be matched in such a way to get a solution which is uniformly valid in the whole flow domain. The matching is carried out using the matching principle of Van Dyke [1], that is, the limit of the outer solution when approaching the corner point must be equal to the limit of the inner solution when approaching infinity. In principle higher approximations are possible by retaining further terms in the asymptotic expansions. But, in practice achieving higher order approximations are difficult due to the amount of calculations involved. In this paper we are satisfied by keeping the leading order terms of the inner and outer solutions.
There are various numerical and experimental studies on dam break flows. One of the closely related study of this problem is the experimental study of Stansby et al. [2] which observes the singularity and shows that horizontal jet occurs at the corner point for classical dam break problem at small times. They also studied this problem for wet-bed case and observed that the mushroom-like jet occurs at the triple point both experimentally and numerically.
Numerical treatment of the dam break flow for small times is rather difficult due to the logarithmic singularity at the bottom corner point. Stansby et al [2] had to artificially remove the singularity by smoothing the corner before applying the numerical method. The numerical treatment of the original dry bed dam break problem is only possible with the help of an asymptotic analysis; near the corner point inner region problem should be solved and the size of the inner region should be provided.
Another related problem of impulsive acceleration of a wavemaker is thought to model the initial stages of the motion of a dam under earthquake loading. The linearized theory of this problem has been studied by Chwang [3]. He solves, to first order, the initial value problem in a small-time expansion and shows the free-surface to be singular at the intersection point. The remarkable point of these flow problems is the jet which occurs at the intersection point with the jet strength being dependent on the angle between the body surface and the free surface of the liquid and the direction of the body motion. This problem is studied for vertical plate in [3 - 5] and for inclined plate [6].
Dam break flows due to gravity is also studied by Pohle [7] and Stoker [8] using the Lagrangian description. However the Lagrangian variables are not suitable for the dam break problems (see Korobkin & Yilmaz [9]). Recently, King & Needham [5], solved a relevant problem of a uniformly accelerating plate into a block of fluid by using Eulerian variables.
They used the leading order outer solution to derive the size of the inner region and then employed an integral transform to obtain the inner solution. Korobkin & Yilmaz [9], used complex analytic function theory to solve the dam break flow problem. They needed the second order outer solution to derive the inner region dimensions. In this paper, we aim to use the Eulerian description and the Mellin Transformation to solve the dam-break problem. Regarding the inner region solution, the free surface shape obtained in this paper is almost the same as that of Korobkin & Yilmaz [9] but the employed method is different. The main idea of the paper is to compare the two methods in terms of methodology and the numerical results. The present method uses the Mellin transform and does not require the second order outer solution to derive the inner one, whereas the second order outer solution is needed for the method by Korobkin and Yilmaz [9].
The second order outer solution is needed to correctly capture the behaviour of the fluid flow near the top corner point and close to the singular corner point at the bottom. Special care is required to calculate the form of the free surface near the top corner point where horizontal and vertical free surfaces meet. In [9] it was theoretically shown that the free surfaces at the top corner point meet up at the order of e2 using the second order outer solution where e is a small parameter introduced to indicate the initial stages of the flow. Immediately after the disappearance of the dam, a part of horizontal and vertical free surfaces is expected to spill over to the region x < 0 (see Fig. 1). Hence to calculate the free surface shape it is necessary to use the Lagrangian variables to follow the fluid particles outside the initial fluid region x > 0, —H < у < 0. The second order outer solution could be found by the Fourier series method just as is done with the first order solution. However the second order problem is far more complicated than the first order problem and it is convenient to use the the domain decomposition method. The main idea of this method is to divide the whole fluid domain into suitable sub domains where solutions can be written as infinite series involving unknown coefficients, and then equate the truncated series at collocation points in the intersection of sub domains to derive the unknowns of the problem. To the authors best knowledge second order outer solution for the whole fluid domain has not been obtained before for the classical dam break problem. The domain decomposition method is used successfully by Needham et al [10] to find the first order outer solution in the problem of an inclined plate accelerating into a body of fluid.
Domain Decomposition methods suggest a convenient way to solve the complicated two and three dimensional nonlinear problems numerically by the concept of domain splitting instead of using arduous finite element approximations for the whole domain [11]. As an example Glowinski et al solved the numerical simulation of a transonic flow by the Schwarz-alternating method with the overlapped regions. They showed the efficiency and the stability of these methods by applying them to the several Poisson problems. Cai [12] is also concentrated on one special group of these Domain Decomposition methods using overlapping subdomains and using the software Diffpack. This study states that the convergence of the solution on the internal boundaries ensures the convergence of the solution in the entire solution domain.
A detailed analysis of the Domain Decomposition methods, and of the 'Schwarz method for overlapping domains' which is similar to the one adopted here, is given in Quarteroni and Valli [13], where the mathematical foundations of the different approaches is provided.
In section 1 we formulate the dam break problem using dimensionless variables, and derive the leading order solution. Next we derive the leading order asymptotics in section 2. Also, we deal with discontinuity at the upper corner point by using Lagrangian variables
in section 3. We derive the second order outer solution using the domain decomposition method in section 4. Finally we draw some conclusions in section 5.
1. Formulation of the problem
The Euler equations of fluid motion are used in the fluid domain together with kinematic and dynamic free-surface conditions at the free surfaces and the slip boundary condition at the bottom to model the dynamics of this system. Dimensionless unprimed variables are introduced as follows,
x' = Hx, y' = Hy, y' = Hy, e = Щ,
t'
Н-г,
9
p' = ррдН.
A mathematical statement of the problem can now be written in non-dimensional form as
du dv
dx + dy 0 du du du dp
dt + U dx + V dy dx
dv dv dv dp
dt + U dx + V dy dy 1,
dy dy
v = m + u---, p
dx
dC
u = m + ^, p
dy
v(x, --- 1,t)
1 < У < У(x,t), İ(y,t) < X < ^
0 on у = y(x, t),
0 on ж = £(y,t), 0,
У(x, 0) = C(У, 0) = ° u{x, y, 0) = v(x, y, 0) = °
as X —— то, u,v — 0 and p — — y,
(1)
(2)
(3)
(4)
(5)
(6)
where equations (1) are Euler equations, (2) and (3) are kinematic and dynamic boundary conditions at the horizontal and vertical free surfaces, (4) is the slip boundary condition at the bottom, (5) is the initial conditions which state that the fluid is at rest initially and (6) is the radiation condition at то.
A small-time solution to (1)-(6) may be sought by posing the power series expansions of the unknown variables, horizontal and vertical components of velocity, horizontal and vertical free surface shapes and pressure, in time,
и = u0(x,y) + tu1(x,y) + 0(t2),
V = Vo(x,y)+ tvi(x,y) + 0(t2),
у = Уо(х)+ try (ж) + t2y^(x) + 0(f), (7)
£ = Co(y) + t$,1(y) + t2^2(y) + 0(t3), p = Po(x,y)+ tpi(x,y) + 0(t2)
as t — 0 x = 0(1), where x = (x,y). We find from the initial conditions that u0 = v0 = 0, ho = hi = 0 and ^o = 6 = 0.
1.1. The leading order solution
By substituting the expansions (7) in (1)-(6) and using the Taylor series expansions of unknown functions about у = 0 for the vertical free surface and x = 0 for the horizontal free surface in the boundary conditions, we fix the domain as the semi-infinite horizontal strip, —1 < у < 0,0 < X < <x> and find the following boundary value problem at the leading order as
Ui
dui
dx dpo dx ’
dvi
dy
Vi =
0,
dpo
dy
1,
1 < У < 0, 0 < X < <x>,
2mx,
1
6 = 2uı(0,y),
Vi(x, —1)=0, У2
Po(x, 0) = 0, po(0,y) = 0 as X —— то, ui,vi — 0 and p0 — —y
(8)
(9)
(10)
(11)
Solution to the problem (8)-(11) is quite straightforward and found by the Fourier Series method as
Po(xyy)
У + ]C
n=0
8(—1)ra
(2n + 1)%2 sin
(2n + 1)2 У
e-{2n+l) X,
Ї2(У)
2 (—1)rasm ^(2n + 1)72y
ж^^ 2n + 1
n=0
(12)
(13)
1.2. Analysis of the singularity near the bottom point (0, -1)
1 ( n \
The sum of the series in (13) is Ç2 = — log ( tan —(1 + y)) which exhibits a singularity in the
ж \ 4 /
free-surface elevation as у — —1. It is seen that the infinite sum in (12) converges at the corner point implying that the pressure is well behaved near the corner point, (0, —1) and that the horizontal fluid velocity and the free surface shape are singular at the corner point,
ui(0,y) = -bg
ж
(tan4(1+
The singularity at the point (0, —1) suggests that the expansions (7) are outer expansions to this problem. An inner solution is required in the neighbourhood of the point (0, —1), as t — 0. First, the behaviour of p0 is investigated as (x2 + y2) 2 — 0. The origin is translated to the bottom point by the transformation £ = x and rq = у + 1 and then the standard polar coordinates (p,9) are introduced by £ = p cos в and q = p sin в (see Fig. 2). Using (12) and (13), as p — 0, we have the limiting values of the outer solution
Po
P log P
(
2
— cos в
* )
+ p
sin в+1 —log(4)
2
+—9 sin 9 ж
+ o(p),
б
1 П 1 , 9л
-log +—logrı + 0(rı ) 'K 4 'к
as p = p ^ 0.
(14)
This indicates that p0 is analytic and £2 is singular at the corner point; p0 = 0(plogp) and £2 = О (log p) as p ^ 0. The fluid velocities in the corner region are calculated by the dp0 dpo
equations u\ = — and v\ = — — 1, from (8). Then we find that the horizontal velocity
d£ dp
has logarithmic singularity at the corner point, u\ = 0(logp) and v\ = 0(1) as p ^ 0. In the next section, an inner region solution is developed near the singular point.
2. Leading order asymptotic (inner) solution near the bottom point (0, -1)
The method of this section follows closely that of King & Needham [5]. In order to formulate an inner solution to this problem when £,p = o(1) as t ^ 0, the magnitude of terms in the velocity components in equations (1) are examined. The local analysis at the end of section 1 shows that as p = (£2 + p2) 2 ^ 0;
и = 0(t log p), V = O(t),
p = 0(plog p), Ç = 0(t2 log p).
Thus a typical term kept in (1) is vt = 0(1) whereas a typical neglected term, which represents fluid inertia, is uvx = 0((t2 log p)/p). These two terms are of equal magnitude when t2 log p = O(p). If this is solved iteratively, we find that p = 0(—t2 log t). In this region inertial terms are important and v = O(t), и = 0(t log t), p = 0(t2 log 2 t) and £ = 0(t2 log t). These estimates suggest the following inner variables,
£ _ = _ V
t2 log d T_ t2 log d
и = и, V = V, p = p,
which gives the boundary value problem in the inner region,
du
~dt
dv
~dt
du dv d- + dfj O’
(}+iıOgî)
du
Ж
1
I (- +--1--^ ^
\t t log t)
f2 log t dq '2 1 N dv f 2
t + t log t) ' \ t
t log t J dr] du 1 dp
t2 log t d-’ f 2 1 \ dv
\t + t log t)
1
dv
-v— =
t log tj dq 1 dp
- 1’
t2 log t dr] t2 log tdr]
V
1 du
-------?/ —=
f2 log t d-
1 dv
------V —=
t2 log t d£
(15)
in the domain 0 < q < ж, £ > £(q,t)/(—t2 log t) and subject to free-surface conditions on
Î = C(V’t)/(-t2 log t),
и
d£ /2 1 \d£ ( 1 \ _di
df - \t + t log tjdfj V-t2 log УЖ and P
0.
(16)
Matching conditions (14) are also written in terms of the inner variables and are applied as p = (.p + rj2) 2 ^ ж. So we get the limiting behaviour of the inner solution as p ^ ж, which suggest the following inner region form of the expansions
p = t2 (log t)2pi + t2 log tp2 + o(t2 log t), £ = t2 log t^i + t2^2 + o(t2),
U = t log tu1 + tu2 + o(t) ’ V = t log tv1 + tv2 + o(t) ’
(17)
(18)
as t ^ 0 with Ç’h = 0(1). If the expansions (17) and (18) are substituted into the inner region problem (15), (16) and by using the Taylor series expansions of unknowns at the free surface, £ = £(fj’t)/(-t2 log t), at the leading order we obtain
du1 dv1 d£ + dq 0’
-du1 du1 du1 du1 dp1
й1 — — 2г>-дП — й1 öf — 51= öp
-dv1 _dv1 _ dv1 _ du1 dp1
V1 - ^ - 2v-Щ — -1 ğ] — -1 ğ- = w j
in the domain 0 < q < ж, £ > —£1, subject to the free-surface conditions on £ = —^,
_d£1 _ d^1
щ = 2& — 2q~— — V1— and p1 = 0. dq dq
From the limiting behaviour of the inner region solution as p ^ ж, at the leading order, we
have,
4 v 2 4
P1 ~ -Ç’ 6 ~ - ’ U1 ~ - ’
ж ж Ж
v1 ~ 0 as (<p + q2)2 ^ ж.
These matching conditions must be the solution or a part of the solution of the leading order problem. By this fact, the exact solution to this problem is
4
Щ = -, Vi = 0,
ж
pi = 4(v + 2), e
ж \ ж)
2
n = -, ж
which represents a mass of fluid moving to the left at the bottom, similar to the findings of [4] where a block of fluid is rising.
Similarly, using the next order terms in the perturbation process, the second order problem is obtained,
дщ dv2 _ 0 5^ + dp ,
4 -du2 du2 4 du2 dp2
- + U2 - 2Ç—2 - 2^—2 - -, ж oÇ dp ж dç dç
_ v^^2 n_dV2 4 dV2 dp2
V2 - 2^Ж - 2^W - _ ~dp - ,
(19)
v 2 v 2
to be solved in the fixed domain 0 < p < ж, £ >-------. The free-surface conditions on £ _------
ж ж
are
4 2 _д&
P2 _ —'^, U2 _ — + 2^2 — 2TW
Ж Ж dp
and the matching conditions, which are obtained from (14), at the second order are,
(20)
P2 ~ ^ {p - 1)C + І log P - ^^, 6 ~ 4(h + log V), (21)
(* - *) 4
22
U2 ~ -p - 1 + log V, V ~-9,
ж ж
(22)
“ (V + v2) 2 - ж ^ „ _ log(- bg t) + Us | 4 |. НЄ.Є M) are the u,lal pobr
coordinates with respect to the Cartesian coordinates (£,??).
Since horizontal velocity u2 and pressure are unbounded as p ^ ж, the boundary value problem (19)-(22) is not easy to solve. This problem is solved in the next subsection.
(0
2.1. Reformulation and the Mellin transform of the boundary value problem (19)
It would be more convenient at this stage to continue with velocity potential p. Since the flow is irrotational at the begining, by the Cauchy — Lagrange theorem, the flow is irrotational
4 dp dp
at all times. Thus a velocity potential p is introduced by u2 +— _ —, v2 _ — and
ж dx dy
- 2
the coordinates are shifted by xi = £ +— and yi = By using this velocity potential,
Ж
the quarter plane problem (19)-(22) becomes
Ay = 0 in 0 < yi < то, x\ > 0, dp 4
Po + 3p - 2yx— = -£ - yi on х-ı = 0, oyi ж
!■
dp 6 df
V— = —+ 2£ - 2yi— on Xi = 0, dx1 ж dy1
dp
— = 0 on yi = 0, oyi
Ф ~ ~\X1 log(xi + vD2 - Vi tan i ( ^ ) + (p +
i ^ + (^ +1)жі j +0 (й + у2)2)
І (logШ + h) + o(1), as n = (xj + УІ)2 ^ TO.
Ж
(23)
An integral transform can not be applied to (23) for ф is unbounded as r\ ^ to. The standard polar coordinates (r\, 0\) relating to the cartesian coordinates (х\,у\) are used and the terms that becomes unbounded as r\ ^ to are subtracted from p and £,
2
p = p +— \г\ cos 9i log Г\ - 0\Г\ sin 0\ + (p + 1)ri cos Q\} ж
І = £ + “(log n + h),
ж
(24)
dp ж
so that p is harmonic with —— (г\, 0) = 0 and the free-surface conditions on 9\ = — can be
UU\ 2
written from (23) by using (24) as
dp
po + 3p - 2ri — dr\
1 dp
r\ddi
-(f + -(log П + ,
^ [ ж J
di
-2« + 2r‘an
It is easily seen that І = o(1) and ф = o(ri) as r\ ^ to. The improper integral of the transform may still be divergent since the behaviour of p at to is o(ri). Thus a coordinate expansion of the form p = A log r\+B + 0(1/r\) for large r\ is suggested, where A = 4/(3^2),
4 8 p0 ~
В = —-p +-----------and І = 0(1/r2). Now using the coordinate expansion, the potential
3ж2 9ж2 3
is redefined as
Ф = Ф +
4
6^2
log(1 + r2)
B,
i= C
Thus we get Ф = 0(іДі) as ri ^ ж and the following boundary value problem,
where
4
'K
Аф = — ——- A log(1 + r-) in 0 < ri < <x>, 0 < Qi < -, 6жг 2
дф
— = 0 on 0i = 0, ов1
дф 4 ж
3Ф — 2Гі— —С + f (ri) on 0i = -,
or1 ж 2
1 дф д( ж
~йГ= — 2^ + 2ri^ on 9i = ö, ri dui ori 2
4 Ґ 3 o 2r?
f * (ri) = 3log ri— 2log(1 +r-) + T+7f
2
)
Now the Mellin transform,
(25)
$Mi)
ri ^(ri,di)dri,
Z (s)
ri i((n)dri,
can be applied to the boundary value problem (25). Since ф = 0(1), ( = О (log ri) as ri ^ 0 and Ф = 0(1/ri), ( = 0(1 /r-) as ri ^ ж we expect Ф to exist and be analytic in the strip 0 < Re(s) < 1 of the complex s-plane and Z will similarly be analytic in 0 < Re(s) < 2. Taking the Mellin transform of the first equation in (25), we have
{.дв\ + s* }Ф(^
(MO = —і
2s
3^ sin(s^/2)
The general solution to this ordinary differential equation is
Ф(з, 0i) = a(s) sin(s0i) + b(s) cos(s0i) —
2
3жз sin(sж/2)
(26)
дф
Similarly taking the Mellin transform of the boundary condition —— = 0 on 0i = 0 gives the
uui
дФ
transformed boundary condition = 0 on 0i = 0. Application of this boundary condition
to the general solution (26) gives a(s) = 0. By transforming the free-surface conditions on ж
di = —, we obtain b(s) as a solution of a difference equation,
b(s)
b(s)(~1)s Г(а)
^sT(s + 5/2)T(s + 2),
where b(s) is a solution of b(s)/b(s — 1) = 1. This gives the Mellin transform of the free-surface elevation by the Mellin transformed free-surface condition (third equation in (25)),
Z(s) = —
b(s)(—1)s 1 cos(s^/2) 2s(s + 1)^s_1 r(s + 3/2)
In order to ensure the convergence of the inversion integral, we require v = Im(s) > —3 and b(s) must have the following behaviour,
0(1), (7 —— + ТО,
0(е2па), а — — то,
b(u + га) = I
where s = V + га. A function of period 1 which has this property is
— c
b(s) = . . -----,
(—1)s sin ns
where c is an arbitrary constant. With this choice of b(p) and by using the following Mellin inversion formulas of Z(p) and Ф(р,в1),
т+гоо т+гоо
СЫ = 2— [ r-sZ(s)ds, ф(гі,ві) = 2~ f r■-sФ(s,вı)ds,
we have,
С (П)
сп 1
т+іоо
(1/пгх )6
4 2пг J s(1 + s) sin(s^/2)r(s + 3/2)
ds
for 0 < т < 2, and
Ф(п,ді) = —
1
2пі
А+гоо
іП
С cos s6ı
А — ІОО
hs sin(^s)s(1 + s)r(s + 5/2) 3^ sin(s^/2)
!
(27)
г- Sds
T — LOO
т—гоо
T—LOO
2
for 0 < X < 1.
2.2. Shape of the free surface and the evaluation of the line integral (27)
Of particular interest now is the form of free-surface that this integral solution represents. The line integral (27) may be turned into a contour integral in Re(s) > 0 by noting that the
-в n n
integrand decays on the semicircle s = pe , — — < в < — and a simple application of residue theorem gives
C (o ) = — 2 ^
(—1)n(1/nn)
2 n
2^2n(2n + 1)Г(2гс + 3/2)
(28)
The ratio test reveals that this series is convergent for all г\ = 0. In fact for large г\ the series is asymptotic and ( = 0(1/^) from (28). However at r\ = 0 the series in (28) diverges
and a rather different approach to the evaluation of ( is needed for small r1. Since the line integral can not be made into a contour integral by addition of a semicircle in the left-s-plane (due to the growth in the gamma function) we consider a rectangular contour surrounding only the poles at s = -1, 0 of the integral (27). Using that rectangular contour, we obtain,
„ 4 CK
с ы = т
{^r(3/2) (log (яті)
1 Г(3/2) + n2n
Г(3/2) + 4
+ o(ri),
(29)
as r1 ^ 0.
Recalling that the physical free-surface elevation is £2 = ( + (1/^)(log r1 + y), the logarithmic term can be eliminated by choosing c as 1Д/7. Thus the form of free surface ((r1) for small r1 in (29) is obtained
U. . . . . n\ Г(3/2) ^2f1
6 = -{ log(-log i) + log(i) - 1 - г(3/2)' 4
+
)
+ o(rı),
and the form of free surface ((r1) for large r1 is
(30)
6
1
Ж
{
log Г1 + log(- log t) + log
E
n=1
(-1)га(1/^Г1)2га
2n(2n + 1)Г(2гс + 3/2)'
(31)
Now, we will show that the outer and inner solutions match satisfactorily using numerical calculations. The correction to the leading order free-surface elevation, £2(rj) can be written as
6(h)
1
-(l°g j + h) +
'K
т+г^
Аг __________________(1Ah6____________d
4 2жі J s(1 + s) sin(s^/2)r(s + 3/2)
т—i<x>
(32)
with 0 < т < 2. The integral appearing in (32) is of 0(1/?j2) as j ^ ж so that £2(j) matches the outer solution as j ^ ж satisfactorily with a relative error of 10-4 as у varies between -0.9 and -0.995 as we see from the Table 1. The series in (31) are calculated for N = 70, which gives a maximum relative error of 10-9.
The light gray region in Table 1 shows the matching between the outer solution and the inner solution for large r and the dark gray region shows the matching between the two forms of the inner solutions.
To compare the expressions for free surface functions (30) with that of Korobkin & Yilmaz [9], it is convenient to rewrite the expression (30) in dimensional variables as,
/
X =
+
gt'2 gt'2 gt'2 ж
----log U7~ +----log 7 -
ж H ж 4
gi
'2
log (- log(f' VgjH))
gt
'2
Ж
Ь + ГШ) +
V + Г(3/2}) +
ж у' + Н
4ю77ўі).
(33)
ж
Table 1. Comparison of the inner solutions (I.S.) of £(y, t) near intersection point for t = 0.04
У Outer Solution I.S. for Large r I.S. for Small r
1 ---0.00034340 ---0.00038318 0.14248134
p
1 ---0.00057253 ---0.00058969 0.09368181
p
1 ---0.00093849 ---0.00094270 0.04488228
p
bo
У = ---0.9 ---0.00129467 ---0.00129571 0.020482515
у = ---0.92 ---0.00140869 ---0.00140936 0.01560256
у = ---0.94 ---0.00155550 ---0.00155586 0.01072260
у = ---0.96 ---0.00176222 ---0.00176234 0.00584265
у = ---0.98 ---0.00211536 ---0.00211525 0.00096270
у = ---0.99 ---0.00246841 ---0.00246781 ---0.00147727
у = ---0.995 ---0.00282143 ---0.00281901 ---0.00269726
у = ---0.999 ---0.00364111 ---0.00358329 ---0.00367325
у = ---0.9992 ---0.00375476 ---0.00366690 ---0.00372205
у = ---0.9994 ---0.00390127 ---0.00375409 ---0.00377085
у = ---0.9996 ---0.00410778 ---0.00382576 ---0.00381965
у = ---0.9998 ---0.00446079 ---0.00386599 ---0.00386845
у = ---0.9999 ---0.00481381 ---0.00389300 ---0.00389285
у = ---0.99992 ---0.00492746 ---0.00389740 ---0.00389773
у = ---0.99994 ---0.00507397 ---0.00390249 ---0.00390261
у = ---0.99996 ---0.00528047 ---0.00385408 ---0.00390749
у = ---0.99998 ---0.00563349 ---2.8631 X 1017 ---0.00391237
у = ---0.99999 ---0.00598651 ---0.00391481
The corresponding expression of Korobkin & Yilmaz [9] is
' gt'2 gt'2 gt' 2 ж gt'2 (1 a\
X = log + log - - log +
ж H ж 4 ж 3 с
gf
'2
(У + Н)2 e
2ж
б 1 gt'2 а ж '2
+ (У + Н) arctan ' - gt S '2
а 3жу + Не 2 ад t2
+
1
9^2
+
СУ + Н \ ад t'2
where
СЮ
S (и) = 2 J ((3 — т 2)e-T /6 — 3cos cos (ur2 /ж) —3, 0
(34)
a = —6 log a, e = gT2/H,
and T is a suitable time scale. The first two terms in the expression (33) match exactly with those of (34) and the third term of (33) is almost the same as that of (34). However it is difficult to compare the last two terms.
Fig. 3. The comparison of the shape of the free surface using the two methods for dimensional variables t = 0.0143 s, H = 1 m, g = 9.81 m/s2
The tangent of the angle between the free surface and the bed is calculated by differentiating the equations (33) and (34). In the leading order both equations (33) and (34) provide the same formula,
dy_
dx'
2
log
ж
(4).
(35)
at the intersection point. Next we compare the expressions numerically and visually in a graph.
In Fig. 3 the shape of the free surface in the inner region is plotted with the dimensional variables using the present method and the method of Korobkin & Yilmaz [9]. Note that the present method which is based on the Mellin transform, is different from that of Korobkin & Yilmaz [9] which makes use of complex analytic function theory. The shapes of free surfaces in Fig. 3 are quite close to each other. Maximum relative discrepancy between the results is about 0.01.
One can observe from the shape of the physical free surface in dimensional variables, shown in Fig. 3, that a jet is formed near the intersection point. Note that very close to the bottom the free surface shape becomes almost perpendicular to the ж-axis. In Fig. 3 the angle between the free surface and the bottom is about 75 deg and it decreases as time increases (see equation (35)). There is also the problem of discontinuity at the point where the horizontal free surface meets the vertical free surface, (0, 0), which we address in the next section.
3. The shape of the free surfaces at the upper corner point
To describe the shape of the free surface at the upper corner point (0,0), Lagrangian variables should be used, for we should be able to follow the fluid particles outside the domain, x > 0. The velocities of the fluid can be written by differentiating the leading order pressure (12),
Fig. 4. The shapes of the free surfaces in dimensionless variables at the leading order for times t = 0.2, t = 0.4 and t = 0.6
dx 8(- 1)n / / n\ ( n \ _(2 K
и = — = -t ^ (2n\ /)2^ - (2n + 1)2 sin (2n + l)-y e ( n+)21, t> 0, (36)
«=S=- (2;;+1)2,4(2«+i)acod(2«+i)av-(2n+i)ғ t>o.
n=0
(2n + l)%2
n=0
8(-l)n
((2n + l)0 cos ( (2n + l)2tfj
The sum of the infinite series in (36) is found as
dx
dt
2t
log
(1 - е“7X)2 + ^2 cos ^ye~7x/2 j
'K
1 + e 7X + 2 sin ye 7x/2
)
2
t> 0,
(37)
dy 2t
dt n
arctan
^ —7X/2^ 2 cos ye x
2y
1- e
V
t > 0.
/
2
The solution to the coupled nonlinear differential equations (37) is carried out by adaptive step size Runge — Kutta numerical routine with an initial condition imposed at t = 0. The shape of the free surfaces near the upper corner point (0,0) is seen in Fig. 4, for different times, with dotted lines denoting the initial shape at t = 0.
4. Second-order outer solution
In order to derive the second-order solution, the domain decomposition method is used [11-13]. In this method, the flow domain is divided into a number of overlapping subdomains. General solutions of the original problem are obtained with some undetermined coefficients in each of the sub-domains. Then the solutions are matched using some collocation points along the boundaries of the sub domains to determine the unknown coefficients. The sizes of the sub domains, the number of the collocation points on the boundaries and the number of terms retained in the solutions of the problem in each sub-domain are the important parameters of the domain decomposition method.
We shall first solve the leading order problem of the dam-break flow by the domain decomposition method and compare the obtained results with the solution (12) by the Fourier series method to determine the optimum parameters for the method. In this method, we divide the flow domain into three regions: Region 1 is the quarter disk of radius r1 about the upper corner point, Region 2 is the quarter disk of radius r2 about the bottom corner point and Region 3 is the semi-infinite strip x > x3 (see Fig. 5). The idea is to find general solutions of the boundary problem (8)-(11) in these three regions. Note that each of the regions include only two of the three boundaries of the total region. Then the solutions are matched on the boundaries of sub regions and the corresponding coefficients of the expansions in each region are determined.
Let r and 9 be the polar coordinates about the origin, x = 0, у = 0, and r and 9 be the polar coordinates centred at x = 0, у = -1. It is convenient to introduce the shifted coordinates denoted by tilde, ж = x, у = у + 1.
A general solution of (8), (10) in Region 1 (R1) is
Po = C„r2n sin(2nd), — < 9 < 0, r < ri.
n=1
A general solution of (8)-(10) in Region 2 (R2) reads,
Po = -У + -n
I1 — log (y)
— J J X — X log r +
+ ^ C^1 r2n 1 cos((2n — 1)r),
n=1
2
'K
0 < 9 < -,
< < 2’
r < t 2,
Fig. 5. Sketch of the decomposition adopte
and, a general solution in Region 3 (R3), which satisfies (8)-(11) except the condition at
X = 0, is
Po
(1 - У) + ^ CnH
cos
n= 1
|^2(2n - i)^j exp (^-2(2n - i)5j ,
X > x3.
The decomposition is applied with r2 = r1, x3 = r1/2 and r1 > 1/2. The coefficients in the expansions are obtained by a collocation method on a number of points which is larger than the number of the unknown coefficients. So that an over-determined system of equations is obtained and solved by the method of Singular Value Decomposition.
Let Ni denote the number of terms in the expansion considered for region Ri, і = 1, 2, 3. The number of collocation points on the boundary of the regions Ri and Rj is chosen to be NcNi + NcNj where Nc is an integer larger than unity. Due to the peculiarity of the solution in region 2, the number of terms for this region is chosen to be two times larger than the other two regions. So, there are three parameters to be selected: the number of terms in the expansions Ni, the free parameter Nc and the radius r1. We get the best solution for the leading order using the parameters N1 = N3 = 8, N2 = 16, Nc = 8 and r1 = 0.8. In order to validate the procedure, the domain decomposition solution is compared with the analytical solution (Table 2).
In Table 2 the pressure at the line segment x = 0.4, -0.98 < у < -0.1 are compared and it is found that the maximum relative error is 2.036 • 10-6. In Table 3, we aim to investigate the effect of the free parameter Nc on the coefficients. It is seen that increasing Nc from 8 to 10 or 12, which has the effect of increasing the number of collocation points, does not change the coefficients significantly.
It is seen from Tables 2 and 3 that the domain decomposition method works well for the leading order problem. Now this method is applied to the second-order problem.
Table 2. Comparison of the leading order pressure po(0A,y) at the distance x = 0.4 from the vertical free surface given by the analytical solution (12) (second column, A.S.) with the solution by the domain decomposition method (third column, D.D.S.)
У A.S. D.D.S. Rel. Err.
у = -0.1 0.03772292 0.03772284 2.036 • 10-6
У = -0.2 0.07618666 0.07618652 1.838 • 10-6
1
р 0.11619080 0.11619080 2.012 • 10-8
оо
1 0.15865982 0.15866034 3.285 • 10-7
р
У = -0.5 0.2047235 0.20472455 5.050 • 10-8
1 0.25581418 0.25581716 1.167 • 10-7
р
1 0.31376662 0.31376923 8.324 • 10-8
р
1
р 0.38085756 0.38085983 5.978 • 10-8
00
У = -0.9 0.45962375 0.45962585 4.558 • 10-8
у = -0.92 0.47698906 0.47698698 4.348 • 10-8
у = -0.94 0.49492189 0.49492395 4.160 • 10-8
у = -0.96 0.51343904 0.51344109 3.988 • 10-8
у = -0.98 0.53254616 0.53254820 3.832 • 10-8
Table 3. Effect of the free parameter Nc on the coefficients
c Nc = 8 Nc = 10 Nc = 12
c{ ---0.589675 ---0.585556 ---0.584781
cl 0.161116 0.161908 0.161601
cl 0.044510 0.044659 0.044744
cl 0.100245 0.100737 0.101034
cl 0.074988 0.075511 0.075808
cl 0.136470 0.138244 0.139307
cl 0.158118 0.160940 0.162633
cl 0.200447 0.200983 0.201161
cl1 0.17168 0.170632 0.169947
cl1 0.152807 0.153092 0.153275
cl1 0.061146 0.061409 0.061579
cl1 0.042178 0.042119 0.042081
cl1 0.033989 0.034421 0.034701
cl1 0.049782 0.049885 0.049947
cl1 0.049615 0.049748 0.049828
cl1 0.023350 0.025080 0.026204
1--1 0.064207 0.064415 0.064544
'“■ч 05
О
1---1 О 0.105125 0.105555 0.105819
1--1 Т-H
О
С[[ 0.049159 0.054885 0.058603
'“-ІСМ 0.148866 0.148929 0.148907
1--1 1---1
О
1---1 СО
1--1 і---1 0.181092 0.184229 0.186180
О
ll
°14 ---0.052022 ---0.030383 ---0.016215
ell 0.508789 0.504884 0.502257
r<II 0.718653 0.732841 0.741539
°16
с111 ---0.414390 ---0.415800 ---0.416708
1--1
1--1 1.525835 1.527664 1.528866
'“-ıCM
О
1--1
1--1 ---1.864034 --- 1.865337
1---1 СО
О 1.866100
с!11 6.345436 6.356601 6.363730
с{и ---13.89871 --- 13.92205 ---13.93554
Г<Ш 49.18452 49.38880 49.51692
°6
с111 ---105.4966 --- 105.5670 ---105.5876
Г<ш 425.9344 427.5115 428.4633
и8
The second-order problem is derived from the boundary value problem
(1)-(6),
△p2 = — 2(u\x + u\>y), —1 < у < 0, 0 < X < ж,
p2(x, 0) = 2'п‘2(х) + Щ(x), P2(0,y) = 2^2(У),
p2,y(%, — 1) = 0, p2 ^ 0 as X ^ ж,
where ^2(ж) and £2(y) are the leading order terms of the horizontal and vertical free surfaces respectively, as given in (9). General solution in Region 1 is written in the form:
P2 = A + (v\ - u\) ( -в + A + A1 log(r)(l>! + 1) + -0^.
2 'к 2 'к 'К 'K
'K
+ CInr2n sin(2nö) - u\, --< в < 0,
n=1
in Region 2,
P2 = 2(u? - V21 + 1) + Y, Cn f2n- cos((2n - 1)0) - «2, 0 < < -,
in Region 3,
n=1
1
P2 = 2(^1 + ^ - uj + uıX + г>1 у)+
+ ^ C™ co^2(2n - 1)^ exp ^2(2n - 1)^ - u\, x > x3.
The velocities of the fluid at the second order can be written by using these second order solutions as
и
V
dxt3
— = tu1 - — (U1 U1,x + г>1 щ,у + p2 ,x), dyt3
= tV1 - — (U1 V1,x + ^1 V1,y + p2 у),
t> 0, t > 0.
As with the first order solution, the solution to the second order problem is carried out by adaptive step size Runge — Kutta numerical routine with an initial condition imposed at t = 0. The second order solution is obtained using the domain decomposition method described above and the comparison of the shapes of the free surfaces near the upper corner point using leading and second order solutions with Lagrangian variables is given in Fig. 6.
Fig. 6. The shapes of the free surfaces in dimen- Fig. 7. The shapes of the free surfaces in dimensionless variables with the comparison of the sionless variables at the second order for differ-leading order and second order for t = 0.2 ent times
Fig. 8. The shapes of the free surfaces at the second order for t
0.2
It is seen that the second order solution makes a larger difference in the vertical free surface than in the horizontal free surface. In Fig. 7, the shapes of the free surfaces near the upper corner point is plotted using both leading and second order solutions for different times. Finally the complete picture of the shapes of the free surfaces using Lagrangian description for the upper part and Eulerian description for the bottom part at the second order can be seen in Fig. 8.
5. Conclusions and remarks
The linear leading-order outer problem of the two-dimensional dam break flow has been solved by the Fourier series method and the Domain decomposition method. Comparison of the leading order solutions of both methods show good agreement. By the analytical solution, it is observed that the flow is log singular at the bottom point which appears as ajet formation starting from this point. The leading-order inner region solution and the correction to the leading order is obtained by using an integral transform method (Mellin transform). From the free surface shape in the inner region it is found that a jet is formed. The comparison of the free surface shapes by the present method and the method of Korobkin & Yilmaz [9] shows that the two methods give quite close results. However due to the particularity of the present method an artificial discontinuity of the free surface appears at a point very close to the bed. The angle between the free surface and the bed predicted by both methods is the same at the leading order.
The second order outer solution is calculated to determine the shape of the free surface at the upper corner point (0,0). Lagrangian variables is needed to follow the fluid particles outside the fluid domain instead of Eulerian variables. Comparison of the shapes of the free surfaces near the upper corner point with leading and second order solutions showed that the second order solution outer makes a larger difference in the vertical free surface than in the horizontal free surface, compared with leading order solution.
The complete picture of the shapes of the free surfaces using Lagrangian description for the upper part and Eulerian description for the bottom part at the second order is obtained.
Acknowledgements. This study was carried out with support from TUBITAK (The Scientific and Technological Research Council of Turkey), Grant No. 111M209. The authors are also grateful to the referees, who suggested several important improvements to the presentation of the results.
References
[1] Van Dyke, M. Perturbation methods in fluid mechanics. New York: Acad. Press, 1964. 229 p.
[2] Stansby, P.K., Chegini, A., Barnes, T.C.D. The initial stages of dam-break flow // J. Fluid Mech. 1998. Vol. 374. P. 407-424.
[3] Chwang, A.T. Nonlinear hydrodynamic pressure on an accelerating plate // Phys. of Fluids. 1983. Vol. 26 P. 383-387.
[4] Peregrine, D.H. Flow due to vertical plate moving in a channel (unpublished note). Department of Mathematics, University of Bristol, U.K. 1972.
[5] King, A.C., Needham, D.J. The initial development of a jet caused by fluid, body and free surface interaction. Part 1. A uniformly accelerating plate // J. Fluid Mech. 1994. Vol. 268. P. 89-101.
[6] Roberts, A.J. Transient free surface flow generated by a moving vertical plate // Quart J. of Mechanics & Appl. Math. 1987. Vol. 40. P. 129-158.
[7] Pohle, F.V. The Lagrangian equations of hydrodynamics: solutions which are analytic functions of the time: PhD Thesis. USA, New York Univ., 1951.
[8] Stoker, J.J. Water waves. New York: Intersci. Publ., 1957. 567 p.
[9] Korobkin, A., Yilmaz, O. The initial stage of dam-break flow // J. Eng. Math. 2009. Vol. 63. P. 293-308.
[10] Needham, D.J., Chamberlain, P.G., Billingham, J. The initial development of a jet caused by fluid, body and free surface interaction. Pt 3. An inclined accelerating plate // Quart J. of Mechanics & Appl. Math. 2008. Vol. 61. P. 581-614.
[11] Glowinski, R., Dinh, Q.V., Periaux, J. Domain decomposition methods for nonlinear problems in fluid dynamics // Comput. Methods in Appl. Mech. and Eng. 1983. Vol. 40. P. 27-109.
[12] Cai, X. Overlapping domain decomposition methods. Advanced topics in computational partial differential equations // Lecture Notes in Comput. Sci. and Eng. 2003. Vol. 33. P. 57-95.
[13] Quarteroni, A., Valli, A. Domain decomposition methods for partial differential equations. New York, UK: Oxford Sci. Publ., 1999. 376 p.
Received 20 December 2018
