Vortex Sheet Formulations and Initial Value Problems:
Analysis and Computing
David Ambrose
August 6, 2014
Goals
Part I:
Describe an efficient numerical method for the solution of the
initial value problem for interfacial flow with surface tension in 2D.
Describe a well-posedness proof for interfacial flow with surface
tension in 2D.
In the water wave case, demonstrate that the zero surface tension
limit can be taken.
Part II:
Do the same for 3D flows.
Also discuss numerical analysis for the 3D numerical method.
Part III:
Describe a new formulation for the traveling wave problem which
allows for overturning waves.
Discuss analysis and computing of traveling waves using this
formulation.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
1 / 67
Part I:
Initial Value Problems in 2D
Includes joint work with Nader Masmoudi.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
2 / 67
Interfacial Flow
We consider two infinitely deep, horizontally periodic fluids,
separated by a sharp interface.
The fluid velocities are given by the irrotational, incompressible
Euler equations:
ρi (vi,t + vi · ∇vi ) = −∇pi ,
div(vi ) = 0,
vi = ∇φi .
The fluids have densities ρ1 and ρ2 .
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
3 / 67
Comments on the Problems in 2D
If the upper density is zero (ρ2 = 0), this is the water wave case.
This problem is well-posed with or without surface tension.
If both densities are nonzero, then the problem is ill-posed without
surface tension (this is related to the Kelvin-Helmholtz
instability). The ill-posedness has been proved by several authors:
Caflisch-Orellana, Duchon-Robert, Lebeau, Kamotski-Lebeau, Wu.
With surface tension, the problem is well-posed regardless of the
densities.
One should only expect to be able to take the zero surface tension
limit in the water wave case.
With surface tension, there is a 3/2-order stiffness constraint from
an explicit timestepping method. Surface tension enters the
problem through the Laplace-Young jump condition for the
pressure: [p] = τ κ, where τ is the surface tension coefficient and κ
is the curvature of the interface.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
4 / 67
Equations of Motion
In either fluid region, the velocity potential is harmonic. Also, the
normal derivative of the potential (i.e., the normal velocity of the
interface) is continuous across the free surface.
This suggest representing the potential with a double layer
potential:
Z
µ(α)(x − x(α), y − y(α)) · n
ˆ (α)
1
dα.
φ(x, y) = −
2
2π
(x − x(α)) + (y − y(α))2
Taking the gradient gives the velocity in the interior (γ = µα ):
Z
1
γ(α)(−(y − y(α)), x − x(α))
(u(x, y), v(x, y)) =
dα.
2π
(x − x(α))2 + (y − y(α))2
We use the Plemelj formulas to take the limit as (x, y) approaches
the interface:
Z
1
γ(α′ )(−(y(α) − y(α′ )), x(α) − x(α′ ))
γ(α) ˆ
PV
dα′ ±
t.
′
2
′
2
2π
(x(α) − x(α )) + (y(α) − y(α ))
2|(xα , yα )|
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
5 / 67
Equations of Motion, continued
In the case of 2D fluids, we can choose to study the horizontally
periodic case: x(α + 2π) = x(α) + 2π,
y(α + 2π) = y(α).
The velocity integral simplifies to an integral over one period.
This is best seen by complexifying: let z(α) = x(α) + iy(α). Then,
the integral is
Z
1
γ(α′ )
dα′ .
2πi
z(α) − z(α′ )
We rewrite this:
Z ∞
∞ Z 2π(j+1)
X
γ(α′ )
γ(α′ )
′
dα
=
dα′
′)
′)
z(α)
−
z(α
z(α)
−
z(α
−∞
j=−∞ 2πj
Z 2π
∞
X
1
dα′
γ(α′ )
=
z(α) − z(α′ ) + 2πj
0
j=−∞
Z
1
1 2π
γ(α′ ) cot
(z(α) − z(α′ )) dα′ .
=
2 0
2
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
6 / 67
An Equation for γ
If we know γ and we know the position of the curve, then we know
the velocity of the curve.
So, we need an expression for γ. Recall that γ is related to the
velocity potentials: γ = µα , and µ is the jump in velocity potential
across the interface.
The velocity potentials satisfy Bernoulli equations.
So, we can get an evolution equation for γ by (a) writing the
Bernoulli equation in either fluid, (b) taking the jump across the
boundary in these to get the evolution equation for µ, and (c)
differentiating this to get the evolution equation for γ.
We will see the result later; the main thing to know is that the
Bernoulli equation contains the pressure, and the pressure jumps
according to the Laplace-Young jump condition, [p] = τ κ. So, the
surface tension appears in the γt equation.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
7 / 67
The Hou-Lowengrub-Shelley Method
Hou, Lowengrub, and Shelley (HLS) introduced a non-stiff
numerical method for 2D interfacial flow with surface tension.
This involves making an arclength parameterization of the free
surface, and computing using (θ, sα ) instead of Cartesian
coordinates (x, y) :
yα
−1
,
s2α = x2α + yα2 .
θ = tan
xα
The free surface has velocity (x, y)t = U n
ˆ + Vˆ
t; from this,
evolution equations for θ and sα can be inferred.
U is determined by physics, but V is chosen to maintain a
favorable parameterization (e.g., arclength) from the equation
sαt = Vα − θα U.
U is decomposed as its most singular part plus a remainder.
With these choices, HLS are able to use a semi-implicit
timestepping scheme and remove the stiffness constraint.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
8 / 67
The θt and sαt Equations
We have the generic evolution (x, y)t = U n
ˆ + Vˆ
t, and the
ˆ
ˆ
geometric identities tα = θα n
ˆ and n
ˆ α = −θα t.
This implies
(xα , yα )t = (Uα + θα V )ˆ
n + (Vα − θα U )ˆ
t.
p
Since sα = x2α + yα2 , we get
sαt =
(xα , yα ) · (xα , yα )t
= (xα , yα )t · ˆ
t = Vα − θα U.
sα
Since θ = tan1 (yα /xα ), we get
θt =
1
1+
2
yα
x2α
·
David Ambrose (Drexel University)
xα yαt − yα xαt
(xα , yα )t · n
ˆ
U α + θα V
=
=
.
2
xα
sα
sα
Vortex Sheet IVPs
August 6, 2014
9 / 67
Equations of Motion, Revisited: Normal Velocity
The normal velocity is dictated by the fluid dynamics. This
implies that U = W · n
ˆ , where W is the Birkhoff-Rott integral.
Using complex notation, denote the interface as z = x + iy. The
Birkhoff-Rott integral is W = (W1 , W2 ), given by
Z
1
1
′
′
(W1 − iW2 )(α) =
PV γ(α ) cot
(z(α) − z(α )) dα′ .
4πi
2
We have an evolution equation for γ :
t)γ)α
θαα ((V − W · ˆ
+
sα
sα
1
ˆ
ˆ
ˆ
− 2At sα Wt · t + 2 γγα − (V − W · t)Wα · t + gyα .
4sα
γt = τ
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
10 / 67
Equations of Motion, Revisited: Tangential Velocity
As we have said, the tangential velocity is chosen to enforce our
preferred parameterization.
We have sαt = Vα − θα U, and we have specified U.
In the periodic case we would like to use a normalized arclength
L
, at all times.
parameterization: sα = 2π
The length, L, evolves according to
Z 2π
θα U dα.
Lt = −
0
So, we should choose V to satisfy
Vα =
Lt
+ θα U = P(θα U ),
2π
where P is the projection which removes the mean of a periodic
function.
The constant of integration does not typically matter.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
11 / 67
Approximation of the Birkhoff-Rott Integral
The Birkhoff-Rott integral can be approximated well with a
Hilbert transform:
Z
1
1
′
′
(W1 − iW2 )(α) =
PV γ(α ) cot
(z(α) − z(α )) dα′
4πi
2
Z
1
γ
γ(α′ )
1
1
′
′
∼
PV
cot
(α − α ) dα = H
.
′
4πi
zα (α )
2
2i
zα
π
n).
L H(γˆ
π
n − Lπ H(γθα )ˆ
t.
L H(γα )ˆ
In real notation, this can be written as W ∼
Differentiating, this becomes Wα ∼
We have made two different approximations – using the first term
of a Taylor expansion, and pulling functions through Hilbert
transforms. The error from both of these approximations can be
controlled, with estimates for the errors in Sobolev spaces.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
12 / 67
Estimates for this Approximation
Define K[z] by
K[z]f (α) =
1
4πi
Z
2π
f (α′ )K(α, α′ ) dα′ ,
0
where
′
K(α, α ) = cot
1
1
1
′
′
(z(α) − z(α )) −
cot
(α − α ) .
2
zα (α′ )
2
Note that this is not a singular integral – the singularities of the
two terms in K cancel.
If z ∈ H s+1 , then K maps H 1 to H s :
kK[z]f kH s ≤ c1 kf kH 1 exp{c2 kzkH s+1 }.
We also use commutator estimates:
k[H, ψ]f kH s ≤ ckf kH s−2 kψkH s .
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
13 / 67
Summary of the Evolution Equations
The evolution equation for θ is θt =
Uα +V θα
.
sα
Since U = W · n
ˆ , the leading-order term here is θt ∼ Wα · n
ˆ . For
this, we use our approximation.
Similarly, in the γ evolution equation, we have Wα · ˆ
t, and we use
our approximation here as well.
We get the following system of evolution equations:
θt =
2π 2
H(γα ) + P,
L2
2πτ
2π 2 γ 2
θαα +
H(θα ) + Q.
L
L2
Here, P and Q are collections of smooth terms which may be
treated routinely in the energy estimates.
γt =
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
14 / 67
HLS Method: Spatial Issues
The method is pseudospectral: some operations are carried out
after taking the FFT, and some operations are carried out in
physical space.
Derivatives and Hilbert transforms are multipliers in Fourier
space, so they are carried out after taking the FFT:
ˆ
H(k)
= −isgn(k),
∂ˆα (k) = ik.
Mutliplication by functions in real space is carried out in real
space.
For integrals, quadrature is done in real space. Any singular
integrals to be computed utilize the alternating-point trapezoid
rule.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
15 / 67
HLS Method: Timestepping
HLS use the fourth-order implicit-explicit timestepping method of
Ascher, Ruuth, and Wetton:
du
= f (u) + νg(u),
dt
1
∆t
f nonlinear, g linear,
25 n+1
4
1
u
− 4un + 3un−1 − un−2 + un−3
12
3
4
= 4f (un ) − 6f (un−1 ) + 4f (un−2 ) − f (un−3 ) + νg(un+1 ).
Notice that the linear part of the evolution equation is evaluated
at step n + 1; this is the implicit part of the method.
In the θ − L formulation, the equations are semilinear; thus,
treating the linear part implicitly removes the strongest stiffness
constraint.
Our u is (θ, γ); solving for un+1 just involves using the FFT and
inverting a 2 × 2 matrix.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
16 / 67
The HLS Findings
The “linearization” of curvature with these choices reduces the
timestepping constraint from 3/2-order to first order.
This allows for investigations of the behavior of solutions in
different parameter regimes. HLS found three types of behavior.
For small surface tension or large initial data, they found
finite-time singularities, with the interface colliding with itself.
For intermediate values of surface tension or intermediate initial
data, they found growth of solutions without singularity formation.
For large surface tension or small initial data, they found behavior
which appeared to repeat in a time-periodic fashion.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
17 / 67
HLS Results: Strong Surface Tension
From T. Hou, J. Lowengrub, and M. Shelley, The long time motion of
vortex sheets with surface tension, Physics of Fluids, vol. 9, 1997, pgs.
1933-1954.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
18 / 67
HLS Results: Intermediate Surface Tension
From T. Hou, J. Lowengrub, and M. Shelley, The long time motion of
vortex sheets with surface tension, Physics of Fluids, vol. 9, 1997, pgs.
1933-1954.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
19 / 67
HLS Results: Weak Surface Tension
From T. Hou, J. Lowengrub, and M. Shelley, The long time motion of
vortex sheets with surface tension, Physics of Fluids, vol. 9, 1997, pgs.
1933-1954.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
20 / 67
Using HLS Ideas for Analysis: 2D
The HLS evolution equations for the vortex sheet with surface
tension are semilinear, and this is important in allowing the
stiffness to be removed.
This is also quite helpful in proving well-posedness of the vortex
sheet with surface tension, by energy methods (Ambrose ’03).
Similarly, Darcy flow can be studied with this formulation
(Ambrose ’04), and the water wave with surface tension can as
well (Ambrose-Masmoudi ’05).
Other authors have since taken up these ideas
(Guo-Hallstrom-Spirn ’07; Cordoba et al in several papers ’09, ’10,
etc...; Tanveer-Ye ’11, ’12).
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
21 / 67
The Well-Posedness Proof (The Energy Method)
1
Introduce a regularized problem.
2
Show solutions exist for the regularized problem, for any positive
value of the regularization parameter.
3
Establish an energy estimate for the regularized problem.
4
Use the energy estimate to establish that the solutions of the
regularized problem exist on a common time interval.
5
Take the limit of the regularized solutions as the regularization
parameter vanishes.
6
Show that the limit solves the original (non-regularized) problem.
7
Use estimates similar to the energy estimate to establish
uniqueness and continuous dependence on the initial data.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
22 / 67
The Approximation Step
There are many ways to form an approximate system.
One goal when forming an approximate system is: You should
know that solutions exist for the approximate system.
My preferred method is to (a) introduce mollifiers into the
evolution equation, and then (b) use the Picard Theorem for
ODEs to get existence of solutions.
Mollifiers are related to the idea of regularization by convolution;
the function ∂α f can be ‘mollified’ by convolving with an
approximate Dirac mass. Call this Jε ∂α f.
For fixed ε > 0, the operator Jε ∂α is a bounded operator.
A Lipschitz estimate for the right-hand sides of the evolution
equations can be performed, giving existence of solutions for a
very short time (this time goes to zero with ε).
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
23 / 67
The Energy Estimate
The key step in the above outline of a well-posedness proof is the
energy estimate which is uniform in the regularization parameter.
We typically find the energy estimate for the unregularized system
first, and then introduce the regularization in a way that does not
disturb the structure of this a priori estimate.
For the 2D vortex sheet with surface tension, we take
E(t) =
k(θ, γ)k2L2
1
+
2
Z
(∂αs θ)2 +
π s−1
(∂ γ)H∂α (∂αs−1 γ)
τL α
π2
+ 2 2 γ 2 (∂αs−1 γ)2 dα.
τ L
With this energy, we are able to show that there exist positive
dE
constants c1 and c2 such that
≤ c1 exp{c2 E}.
dt
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
24 / 67
Main Idea for the Energy Estimate
The energy controls kθkH s and kγkH s−1/2 .
We take a time derivative:
Z
Z
d1
(∂αs θ)2 dα = (∂αs θt )(∂αs θ) dα
dt 2
Z
2π 2
(H∂αs+1 γ)(∂αs θ) dα + l.o.t. (1)
= 2
L
Note that this cannot be controlled in terms of the energy.
We continue:
Z
Z
π
d1 π
(∂αs−1 γ)(H∂αs γ) dα =
(∂αs−1 γt )(H∂αs γ) dα + l.o.t.
dt 2 τ L
τL
Z
Z
2π 2
2π 3
s+1
s
= 2
(∂α θ)(H∂α γ) dα+ 3 γ 2 (H∂αs θ)(H∂αs γ) dα+l.o.t.
L
τL
(2)
There is an important cancellation when adding (1) and (2),
leaving only one term which is not bounded by the energy.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
25 / 67
Main Idea for the Energy Estimate, continued
We continue:
d 1 π2
dt 2 τ 2 L2
Z
γ
2
Z
π2
dα = 2 2 γ 2 (∂αs−1 γt )(∂αs−1 γ) dα+l.o.t.
τ L
Z
2π 3
γ 2 (∂αs+1 θ)(∂αs−1 γ) dα + l.o.t. (3)
=
τ L3
(∂αs−1 γ)2
Adding (1), (2), and (3), there is again a cancellation. All the
remaining terms can be bounded in terms of the energy:
dE
≤ c1 exp{c2 E}.
dt
Solving this differential inequality, we find
E(t) ≤ −
David Ambrose (Drexel University)
ln(exp(−c2 E(0)) − c1 c2 t)
.
c2
Vortex Sheet IVPs
August 6, 2014
26 / 67
The Well-Posedness Result
Theorem
Given θ0 ∈ H s , γ0 ∈ H s−1/2 (subject to certain conditions), there exists
T > 0 such that there is a unique solution (θ, γ) of the initial value
problem for the vortex sheet with surface tension with initial conditions
(θ0 , γ0 ), with θ ∈ C([0, T ]; H s ) and γ ∈ C([0, T ]; H s−1/2 ). The time T
depends on kθ0 kH s , kγ0 kH s , the physical parameters (such as τ, ρ1 , and
ρ2 ), and on a couple of other things. The solution depends continuously
′
′
on the initial data in H s × H s −1/2 for any s′ < s.
The conditions mentioned above: θ0 must be the tangent angle for
a 2π-periodic curve, and the curve must not have self-intersection
(a chord-arc condition is satisfied).
This is for any values of ρ1 , ρ2 (not both zero) and τ > 0.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
27 / 67
Summary of 2D Problem with Surface Tension
Solutions have been shown to exist for the initial value problem for
the vortex sheet with surface tension, with arbitrary densities.
In the case that the upper density is zero, this means that solutions
have been shown to exist for the water wave with surface tension.
The surface tension force was very much relied upon in the proofs.
In particular, this means that the time of existence for solutions
vanishes as the surface tension parameter vanishes.
In the water wave case, we can do better than this.
Other than the water wave case, we cannot do better – the
problems without surface tension are ill-posed, so we cannot
expect to be able to take the zero surface tension limit.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
28 / 67
The Zero Surface Tension Limit
We restrict ourselves now to the water wave case (ρ2 = 0).
In order to take the zero surface tension limit, we need to make an
estimate for solutions with surface tension which is uniform in τ.
To do this, we have found that a change of variables is helpful: we
replace γ with a quantity we call δ, which we named the “modified
tangential velocity.”
We define
δ=
γ
+ W ·ˆ
t − V,
2sα
and we can interpret this one of two ways: (a) δ is like γ plus a
perturbation, or (b) the quantity δ is the difference between the
Lagrangian tangential velocity of a particle on the free surface and
our artifical tangential velocity.
Instead of evolution equations for (θ, γ), we write evolution
equations for (θ, δα ).
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
29 / 67
The Zero Surface Tension Limit, continued
The leading-order parts of the evolution equations are then
θt ∼ H(δα ),
δαt ∼ τ ∂α3 θ + k(α, t)θα .
This k(α, t) is equal to ∇p · n
ˆ < 0, which (as Wu has shown) has a
fixed sign.
We can make estimates uniform in τ. This is similar to what we
would have if we made the estimate when τ = 0.
If τ = 0, we make the following estimate:
d1
dt 2
=
Z
Z
(−k)(∂αs θ)2 + (∂αs−1 δα )(H∂αs δα ) dα
Z
= (−k)(∂αs θt )(∂αs θ) + (∂αs−1 δαt )(H∂αs δα ) dα + l.o.t.
(−k)(H∂αs δα )(∂αs θ) + (k)(∂αs θ)(H∂αs δα ) dα + l.o.t. = l.o.t.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
30 / 67
Comments on long-time existence
For either gravity waves or capillary waves, there are proofs that
small solutions exist for all time.
One cannot expect large solutions to exist for all time.
In 2D or 3D, such results are by Wu, Germain-Masmoudi-Shatah,
Ionescu-Pusateri, Alazard-Delort, and Hunter-Ifrim-Tataru.
An interesting wrinkle is that the case of long-time existence of
gravity-capillary waves is still open.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
31 / 67
Part II:
Initial Value Problems in 3D
Includes joint work with Nader Masmoudi, Michael Siegel, Svetlana
Tlupova, Yang Liu.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
32 / 67
Description of the 3D Problems
Two 3D fluids, one above the other, separated by a sharp interface.
Horizontally, doubly periodic. Vertically, of infinite extent.
In addition to the irrotational Euler problems considered
previously, we will also discuss the case in which the fluid
velocities are given by Darcy’s Law: Vi = −Ki ∇(pi + ρi gz).
Incompressible: ∇ · Vi = 0.
V1 · n
ˆ = V2 · n
ˆ , but there is a jump in the tangential velocity.
p1 − p2 = τ κ
We let the free surface be given by X(α, β, t). We write
Xt = U n
ˆ + V1ˆ
t1 + V2ˆ
t2 .
Putting this together, in the Darcy problem, we have that Xt is
like three derivatives of X, so we can expect a third-order stiffness
constraint from an explicit numerical method.
In the irrotational Euler problems with surface tension, we again
have a 3/2-order stiffness constraint from an explicit method.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
33 / 67
Extending HLS Ideas for Analysis: 3D
3D is, as usual, more difficult. Need to replace the arclength
parameterization and the variable θ.
One benefit of using θ in the 2D case is that reconstruction of the
free surface is straightforward.
In 3D flows, Ambrose-Masmoudi used κ, the mean curvature, as
the main dependent variable for energy estimates.
Also, an isothermal parameterization was used: if the free surface
is X(α, β), then we require
E = Xα · Xα = Xβ · Xβ = G;
F = Xα · Xβ = 0.
With these choices, energy estimates and well-posedness proofs are
again possible (Ambrose-Masmoudi ’07, ’09, Ambrose ’07,
Wang-Zhang-Zhang ’12, Cordoba-Cordoba-Gancedo ’13).
A difficulty is reconstructing X from κ.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
34 / 67
Vortex Sheet with Surface Tension: Evolution Equations
Xt = U n
ˆ + V1ˆ
t1 + V2ˆ
t2 , with U the normal component of the
Birkhoff-Rott integral.
For the vortex sheet with surface tension, recall that µ is the jump
across the free surface of the velocity potentials, and the potentials
satisfy Bernoulli equations.
We thus have an evolution equation for µ :
µt = τ κ +
t2 )µβ
(V1 − W · ˆ
t1 )µα (V2 − W · ˆ
√
√
+
.
E
E
(This is the density-matched case; if the densities are unequal, we
get several more terms.)
We find an evolution equation for κ :
1
1
1
1
1
κt = √
∂α √ Λ √ ∂α µ + ∂β √ Λ √ ∂β µ + l.o.t.
4 E
E
E
E
E
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
35 / 67
Vortex Sheet with Surface Tension: Energy Estimates
We define A by A =
(A1 , A2 )
We define the operator L by
µα
1/2 µβ
1/2
= Λ √ ,Λ √
.
E
E
1
1
L(B 1 , B 2 ) = ∂α √ Λ1/2 B 1 + ∂β √ Λ1/2 B 2 .
E
E
We get the following system:
1
κt = √ LA + l.o.t.,
4 E
At = −τ L∗ κ + l.o.t.
We can then make energy estimates, with leading-order energy
2
k 2
Z 1 1
1 ∗ 1
k ∗
κ dαdβ.
(L √ L) A + τ √ √ LL
E=
4
E
E
E
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
36 / 67
Vortex Sheet with Surface Tension: Existence Proof
The point of the above: good estimates are available for κ.
For the existence argument, we make an iterative scheme. We
assume that we have in hand (Xn , κn , E n , µn , An ), with some
estimate.
We construct κn+1 and An+1 as the solution of a linear evolution
equation, with coefficients based on the previous iterates.
˜ n+1 as a
Then, Xn+1 needs to be constructed: first, we define X
solution of the evolution equation based on the previous iterates,
and
−∆Xn+1 + Xn+1 = −2Eκn+1
˜ n+1 × X
˜ n+1 )
(X
α
β
˜ n+1 .
+X
˜n ×X
˜ n|
|X
α
β
Estimates for the iterates are established, and so on.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
37 / 67
3D Zero Surface Tension Limit
This is similar to the 2D case.
We focus on making an estimate for water waves with surface
tension which is uniform in the surface tension parameter.
As in the 2D case, we make a change of variables, and are able to
estimate the new quantities, either with or without surface tension.
Existence of water waves with surface tension, plus estimates
uniform in surface tension, allow the zero surface tension limit to
be taken, giving another proof of the existence of water waves
without surface tension.
Everything is harder in the 3D case, but the elements of the proof
are essentially a combination of the well-posedness argument for
the 3D vortex sheet with surface tension and the 2D zero surface
tension limit of water waves.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
38 / 67
3D Numerical Method
For the initial surface, we must find an isothermal
parameterization.
Tangential velocities are chosen to maintain the isothermal
parameterization: the equations Et = Gt and Ft = 0 imply the
following:
V1
V2
U (L − N )
√
,
− √
=
E
E α
E β
V1
V2
2U M
√
+ √
=
,
E
E α
E β
where L, M, N are the second fundamental coefficients of the
surface.
The Birkhoff-Rott integral must be computed; we use Ewald
summation.
We make a Small-Scale Decomposition of the evolution equations,
separating out the most singular terms.
We then use a semi-implicit timestepping method.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
39 / 67
The Birkhoff-Rott Integral
The normal velocity is equal to the normal component of the
Birkhoff-Rott integral:
Z ∞Z ∞
(X − X′ )
1
(µ′α X′β − µ′β X′α ) ×
W=
PV
dα′ dβ ′ .
′ |3
4π
|X
−
X
−∞ −∞
For Darcy flow (the stiffest of the problems), we have
µα = −Bκα − W zα + AEW · ˆ
t1 ,
µβ = −Bκβ − W zβ + AGW · ˆ
t2 .
B measures surface tension, W is proportional to the density
difference, and A is proportional to the viscosity difference.
Notice that if A is nonzero, then we have integral equations to
solve for µα , µβ .
Computing W is much harder for 3D flows than for 2D flows.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
40 / 67
Understanding the Birkhoff-Rott Integral
We approximate the Birkhoff-Rott integral.
We add and subtract in the kernel:
If X ∈ H s+1 , then the operator K[X] (which is the integral
operator with kernel K) maps from H s−3/2 to H s .
This gives a decomposition of the Birkhoff-Rott integral into
operators with straightforward symbols (such as Riesz transforms)
and the operator K[X].
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
41 / 67
The Small-Scale Decomposition and Timestepping
We write
n+1ˆ1n
n+1ˆ2n
Xn+1 − ∆t Usn+1 n
ˆ n + V1s
t + V2s
t
= Xn
n ˆ1n
n ˆ2n
+ ∆t (U n − Usn )ˆ
nn + (V1n − V1s
)t + (V2n − V2s
)t
.
n+1
n+1
Here, Usn+1 , V1s
, and V2s
are the most singular parts of the
velocities. For example:
!!
n+1 n+1
µ
1
µ
β
Usn+1 = −
H1 √α
,
+ H2 √
2
En
En
µn+1 = −B
ˆn
Xn+1
ˆ n + Xn+1
αα · n
ββ · n
2E n
!
− W z n+1 .
We have a linear system for Xn+1 , which we solve with
preconditioned GMRES.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
42 / 67
Results: Removing the Stiffness (Darcy Flow with
Surface Tension)
The above timestepping scheme is first-order in time; higher-order
schemes are available, and we have also implemented a
second-order version.
Using our small-scale decomposition, we are able to effectively
remove the stiffness from the problem.
A fully explicit method would have a third-order stiffness
constraint. We instead face only a first-order stiffness constraint.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
43 / 67
Results: Rayleigh-Taylor with Surface Tension
Irrotational Euler flow.
Heavy fluid above a light fluid, with surface tension present.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
44 / 67
Results: Darcy Flow/Finger Formation
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
45 / 67
Convergence Analysis
We want to prove convergence of (a version of) the numerical
method.
In 2D, there are convergence proofs of the HLS method
(Ceniceros-Hou) and of a boundary integral method for water
waves (Beale-Hou-Lowengrub).
In 3D, there are convergence proofs for boundary integral methods
for water waves (Beale; Hou-Zhang).
We are unaware of such a proof of convergence for a 3D boundary
integral method for interfacial flow with surface tension.
To show convergence, we need to show consistency and stability.
Consistency is fairly straightforward.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
46 / 67
About Stability
So, we want to prove stability for a version of our numerical
scheme.
This means that if we have a continuous surface X, and a
computed, discretized surface Xh , we need to prove an estimate
for the growth of X − Xh .
We actually do this for the semi-discrete system (continuous in
time, spatially discrete).
The required estimates are very much like the energy estimates
that we prove for continuous problems.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
47 / 67
Summary So Far, of Analysis and Computing
For computing, we can develop a small-scale decomposition for
evolution of the free surface.
Analytically, we get good estimates for curvature, but of course we
need the surface as well; we went through a complicated procedure
reconstructing a surface based on the curvature at each iteration.
In the numerical analysis, we have a problem, then: we could make
good estimates for curvature, but we are evolving the surface itself.
Our goal is to find a version of the numerical method for which we
can make the desired estimates; we would like the method to be as
close as possible to what was implemented, but the main goal is to
prove convergence of a boundary integral method for 3D
interfacial flow with surface tension.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
48 / 67
Important Ideas for Stability
To prove stability, we prove energy estimates for the difference of
the continuous solution, X, and the computed solution, Xh .
These energy estimates are similar to the estimates of Ambrose
and Ambrose-Masmoudi for related problems.
Things are more difficult in the discrete setting since relationships
may not hold exactly.
For example, the equation Fh = 0 is not satisfied exactly. Idea:
instead, require ∂t Fh = ǫ∆Fh .
Also, estimates work well for κ, but we need to evolve X. We can
overcome this, but we need to alter the method (as compared to
the method we have previously implemented).
With these two key ideas, plus the framework for estimates from
Ambrose and Ambrose-Masmoudi, we are able to close the
stability estimates.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
49 / 67
Part III:
Periodic Traveling Waves in 2D
Includes joint work with Benjamin Akers, J. Douglas Wright, Walter
Strauss.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
50 / 67
A Traveling Parameterized Curve
The usual traveling wave formulation is to posit that the
dependent variables are functions of x − ct, for some wave speed c.
This is fine for us at small amplitudes, but we should be able to
find large-amplitude solutions as well.
At large amplitudes, our traveling vortex sheets with surface
tension should be able to overturn, i.e., we would like to allow for
waves with multi-valued height.
Instead of studying functions of x − ct, we must make a traveling
wave ansatz for a parameterized curve.
With α corresponding to the normalized arclength
parameterization, we seek solutions for which
(x(α, t), y(α, t))t = (c, 0).
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
51 / 67
The Traveling Wave Ansatz
We have two expressions for the velocity of the curve:
(x, y)t = (c, 0),
(x, y)t = U n
ˆ + Vˆ
t.
We also have an expression for the tangential velocity:
Vα = P(θα U ).
The solution of these equations is
U = −c sin(θ),
V = c cos(θ).
The equation for V is redundant, from the equation for Vα above.
We still need another equation.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
52 / 67
Completing the Traveling Wave Ansatz
Since we have a system of evolution equations (for θ and γ), we
need a pair of equations in our traveling wave ansatz.
In addition to specifying U = −c sin(θ), we also specify
Ut = (−c sin(θ))t .
It’s easy to see that θt = 0, so this becomes Ut = 0.
U is the normal component of the Birkhoff-Rott integral:
iθ
Z 2π
1
e L
γ(α′ ) cot
z(α) − z(α′ )
PV
dα′ .
U = Re
2
2
0
Given the equation U = −c sin(θ), this becomes, upon
differentiating with respect to t,
iθ
Z 2π
1
e L
′
′
′
γt (α ) cot
PV
(z(α) − z(α )) dα .
Ut = Re
2
2
0
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
53 / 67
The Traveling Wave Equations
From the above formula for Ut , we can see that if γt = 0, then
Ut = 0.
We solve the equations U = −c sin(θ) and γt = 0, while using the
normalized arclength parameterization.
More explicitly, these equations are:
iθ
Z 2π
1
e L
′
′
′
γ(α ) cot
z(α) − z(α )
PV
dα = −c sin(θ),
Re
2
2
0
τ θαα
+
sα
(V − W · ˆ
t)γ
sα
= 0.
α
(This is the density-matched case; the general case can also be
treated.)
The quantities L, V, and z are all determined by θ, and W is
determined by θ and γ.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
54 / 67
Local Bifurcation Theory
We use a version of a bifurcation theorem from Zeidler, related to
the Crandall-Rabinowitz “bifurcation from a simple eigenvalue”
theory. The theorem has four hypotheses:
1
2
3
4
ψ : H1 × R → H2 is C 2 ,
∃u0 ∈ H1 such that ∀c ∈ R, ψ(u0 , c) = 0,
∃c0 ∈ R such that the linearization L(c0 ) = ψu (u0 , c0 ) has a
one-dimensional kernel and Fredholm index zero,
If ker(L(c0 )) = span(e) and if ker(L∗ (c0 )) = span(e∗ ), then
he∗ , ψuc (u0 , c0 )eiH2 6= 0.
Then, there exist infinitely many solutions (u∗ , c∗ ) such that
ψ(u∗ , c∗ ) = 0 such that u∗ 6= u0 .
This is not the most general version of such a theorem, but it
sufficient for the present purpose of proving that there exist
periodic traveling vortex sheets with surface tension.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
55 / 67
Applying the Bifurcation Theorem
The solution θ = 0, γ = γ0 is the trivial solution for any c.
If we were to choose H1 and H2 to simply be, for instance, some
Sobolev spaces or C k spaces, then the kernels of the operators
would be too large. We reduce the size of the kernels by imposing
symmetry. Thus, θ will be taken to be odd and γ will be taken so
that ω := γ − γ0 is even.
The linearized operator is
L(c)(θ, ω) =
c
H
τ ∂α2 + 2γ02 ∂α H c∂α
θ
ω
.
Hypotheses (3) and (4) hold by direct calculation, using the
symmetry assumptions.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
56 / 67
The Bifurcation
The speed c must be chosen so that, for some nonzero integer k,
c2 = τ |k| − 2γ02 .
Letting k0 be the least positive integer k for which the right-hand
side is non-negative, the first bifurcation value for c is
c20 = τ |k0 | − 2γ02 .
1.2
1
c
0.8
0.6
0.4
0.2
0
1
2
3
4
5
6
|y|
Numerically, we fix τ and use various values of γ0 .
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
57 / 67
Computed Solutions: γ0 = 0
The computations are performed by finding solutions of ψ(u, c) = 0
with a Newton solver, starting with guesses from the linearization
at small amplitudes, and continuing to large amplitudes.
With γ0 = 0, the computed solutions never overturn.
#%
%
#$
&
%
(
"
$
$
"
%
&
#$
%
#%
!
"
#
$
#
"
!
!
'
David Ambrose (Drexel University)
"
#
$
#
"
!
&
Vortex Sheet IVPs
August 6, 2014
58 / 67
Computed Solutions: γ0 6= 0
For nonzero γ0 , we find that the waves always overturn, and the
bifurcation curve ends in self-intersection of the interface.
!
#&
"%&
#%
"
#"
#%&
#$
#
(
'
$%&
&
$
%
$%&
"
#
$
#%&
"
!
"
#
$
#
"
!
"
!
'
David Ambrose (Drexel University)
"
#
$
#
"
!
(
Vortex Sheet IVPs
August 6, 2014
59 / 67
A Close-up of the Interface
&
#"%
(
#
$"%
$
!"%
!"#
!"$%
!"$
!"!%
!
!"!%
!"$
!"$%
!"#
'
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
60 / 67
Another Traveling Wave
The previous waves had ρ1 = ρ2 ; now, instead, we set ρ2 = 0.
A large-amplitude gravity-capillary water wave:
"
(
#%&
#
$%&
!
"
#
$
#
"
!
'
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
61 / 67
A Close-up of the Interface
1.4
1.3
y
1.2
1.1
1
0.9
0.8
2.35
2.4
2.45
2.5
2.55
2.6
2.65
2.7
2.75
2.8
x
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
62 / 67
Gravity Perturbed Crapper Waves
Water waves with fixed surface tension, and a gravity parameter,
g.
%"#
%
'(h
$"#
$
!"#
!
ï!"#
David Ambrose (Drexel University)
!
!"#
$
$"#
%
&
Vortex Sheet IVPs
August 6, 2014
63 / 67
Our Largest Gravity Perturbed Crapper Wave
%$
&
!
"
(
#
$
ï#
ï"
ï!
ï&
ï%$
ï!
ï"
David Ambrose (Drexel University)
ï#
$
#
"
!
'
Vortex Sheet IVPs
August 6, 2014
64 / 67
Global Bifurcation
In current work, we have proved a global bifurcation theorem for
periodic traveling interfacial capillary-gravity waves on infinite
depth.
The densities can be arbitrary, as long as they are not both zero.
The formulation allows for overturning waves, without using a
conformal map.
This gives hope that the same can be done in 3D.
There are several conditions for how the curve of solutions can
end, including: blowup of solutions, the curve going to infinity in
some sense, self-intersection of solutions, or reconnection to trivial.
The theorem uses an “identity plus compact” global bifurcation
theorem; the compactness is provided by inverting the derivatives
present in the surface tension term in the evolution equations.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
65 / 67
Current and Future Work
Finish writing the proof of convergence of the 3D boundary
integral method.
Possibly extend the numerical analysis to vortex sheets or water
waves with surface tension
Allow for adaptive mesh refinement via overlapping coordinate
patches, in 2D or 3D.
Parallelize the computation of the Birkhoff-Rott integral in 3D.
Apply these tools to other physical problems (including computing,
numerical analysis, and analysis); for example, hydroelastic waves.
Further traveling wave studies, especially the study of overturning
traveling waves in 3D.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
66 / 67
Thanks for your attention.
David Ambrose (Drexel University)
Vortex Sheet IVPs
August 6, 2014
67 / 67
© Copyright 2024 ExpyDoc
