Variational space–time (dis)continuous Galerkin method for nonlinear free surface water waves

(1)

Contents lists available atScienceDirect

Journal

of

Computational

Physics

www.elsevier.com/locate/jcp

Variational

space–time

(dis)continuous

Galerkin

method

for

nonlinear

free

surface

water

waves

E. Gagarina

a

,

∗

,

V.R. Ambati

a

,

J.J.W. van der Vegt

a

,

O. Bokhove

b

,

a

,

∗

a_Mathematics_of_{Computational}_Science_Group,_Dept._of_Applied_Mathematics,_University_of_Twente,_P.O._Box_217,₇₅₀₀_AE,_Enschede, The Netherlands

b_School_of_Mathematics,_University_of_Leeds,_LS2_9JT,_Leeds,_United_Kingdom

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received25October2013

Receivedinrevisedform21May2014 Accepted18June2014

Availableonline26June2014 Keywords:

Nonlinearwaterwaves FiniteelementGalerkinmethod Variationalformulation Symplectictimeintegration Deforminggrids

A new variational finite element method is developed for nonlinear free surface gravity water waves using the potential flow approximation. This method also handles waves generated by a wave maker. Its formulation stems from Miles’ variational principle for water waves together with a finite element discretization that is continuous in space and discontinuous in time. One novel feature of this variational finite element approach is that the free surface evolution is variationally dependent on the mesh deformation vis-à-vis the mesh deformation being geometrically dependent on free surface evolution. Another key feature is the use of a variational (dis)continuous Galerkin finite element discretization in time. Moreover, in the absence of a wave maker, it is shown to be equivalent to the second order symplectic Störmer–Verlet time stepping scheme for the free-surface degrees of freedom. These key features add to the stability of the numerical method. Finally, the resulting numerical scheme is verified against nonlinear analytical solutions with long time simulations and validated against experimental measurements of driven wave solutions in a wave basin of the Maritime Research Institute Netherlands.

1. Introduction

Thehydrodynamicsofwaterwavesinoffshoreregions ofseasandoceansareofsigniﬁcantinterest tonavalarchitects andmarineengineers.Particularlyimportantisthedynamicsofextremewavephenomena,suchasroguewaves,andtheir effectonoffshorestructuresandships.Engineersoftentestprototypemodelsofshipsandoffshorestructuresin experimen-talwavebasinswithwavemakerstomimicrealisticmaritimeenvironments.Forextremewavesitisnontrivialtogenerate waveswithspeciﬁcpropertiesatcertainmeasurementlocationsinawavebasin.Numericalmodelingofwavesistherefore animportantcomplementaryactivity.

A class of water wave problems is described by the Laplace equation for the velocity potential with Neumann and nonlinearfreesurfaceboundaryconditions.ThesefreesurfacegravitywaveequationsareobtainedfromtheEulerequations offluidmotionundertheassumptionsthatthefluidisinviscidandincompressible,andthevelocityfieldirrotational

[18]

.

*

Correspondingauthors.

E-mailaddresses:E.Gagarina@utwente.nl(E. Gagarina),vrambati@gmail.com(V.R. Ambati),j.j.w.vandervegt@utwente.nl(J.J.W. van der Vegt), o.bokhove@leeds.ac.uk(O. Bokhove).

(2)

Free surface gravity waterwave equationscanalso beobtained inasuccinct wayvia Luke’svariationalprinciple

[26]

orfromitsdynamical equivalentpresentedby Miles

[30]

.Inthevariationalprinciple,thecompleteproblemisembedded in asingle functional.In addition,thevariational formulationisassociatedwithconservationofenergy andphase space. Improvedfiniteelementsolutionsofthefreesurfacewaterwaveequationscanbeconstructeddirectlyfromthisvariational formulationratherthanfromaseparateweakformulation.Sincevariationalfiniteelementformulationshaveasmainbenefit thattheyprovideanaturalsettingtoensurediscreteenergyandphasespaceconservation,theycanmorenaturallyleadto stablenumericalschemes.

Variational andﬁnite element formulations for free surface waves usingpotential ﬂow can be found in Kim and Bai

[19] and Kim et al. [20]. Their formulation is, however, restricted to an unbounded domain and involves a coordinate transformation toa referencedomain.Klopmanetal.

[21]

derived avariational Boussinesqmodel fromMiles’variational principle (see also [6,12]). More, standard ﬁniteelement methods for free surface gravity waterwaves can be found in

[27–29,42–45].

Another widely used numericalmethod forfree surface wavesusing potential ﬂow is the boundary integral method, whichstartedwiththeworkofLonguet-HigginsandCokelet

[25]

,andVinjeandBrevig

[40]

.Thismethodhasbeenapplied extensivelytotwodimensionalfreesurface waves,seethesurveysofRomateandZandbergen

[33]

andTsaiandYue

[37]

. Applicationsinthreedimensions(3D) usingboundaryintegralmethodsarewidespread

[3,5,9,10,14,17]

.Acomparison be-tween the workloadofboundary andﬁnite elementmethods ismadeinBettess [4], withﬁnite elementmethods being competitiveinanisotropicdomains.

Morerecently,discontinuousGalerkin(DG)methodsforellipticproblemshaveenabledresearcherstomodelfreesurface waterwaveswithDGmethods(seetheoverviewinArnoldetal.

[2]

).Inparticular,space–timeDGmethods,inwhichtimeis consideredasanadditionalcoordinateandtheproblemissolvedind

+

1 space,withd thespatialdimension,arepromising forfreesurface ﬂowswiththeirdeforming domains.Space(–time)DGﬁniteelementmethodsforlinearfree surfacewave problemsweredevelopedinVan derVegtandTomar

[38]

andfornonlinearflowsinVanderVegtandXu[39].Thelatter approach suffers,however, from instabilities forhigher amplitude waves because the couplingbetween the free surface dynamicsandtheinteriormeshmovementissuboptimal.We haveovercomethisshortcomingbydirectlydiscretizingthe fluiddynamicalvariationalprincipleforpotentialflowwaterwaves.

Numericallysolvingthenonlinearfreesurfacegravitywaveequationsposesanumberofchallenges:

(i) thesolutionofthegoverningequationsdependsonthepositionofthefreesurface,whichisnotknownapriori; (ii) theﬂowdomaincontinuously deformsintime,notonlyduetothefreesurfacemotion,butalsoduetothepresence

ofawavemaker;

(iii) energy drift needs to be avoidedin time inorder tominimize numerical decayofwave amplitude (in the unforced limit);and,

(iv) dispersionerrorsneedtobesmallforaccuratecapturingofcomplexwaveinteractions.

To copewiththesechallenges,we developeda variationalﬁniteelement methodfornonlinearfree surfacegravity water wavesbasedonMiles’variationalprinciple.

In our space-plus-time approach, space is discretized first by considering expansions in terms of continuous basis functions inspace. Theseare substituted directlyintoMiles’ variationalprinciple. Theresulting spacediscrete variational principleiseitherautonomous(i.e.,withnoexplicittimedependence)ornon-autonomous(i.e.,includingtheexplicittime dependenceduetothewavemaker)intime.Subsequently,wedevelopatime(dis)continuousGalerkinapproachbyusing discontinuousbasisfunctionsintimetodiscretizethisspatiallydiscretevariationalprinciple,byfirsteliminatingtheinterior degreesoffreedom. Theresultis thatthefree surface becomes continuouswhile the(surface) velocitypotential remains discontinuous intime. A posteriori,we alsoestablishthetime discretizationoftheinteriordegreesoffreedom. Thefinal resultisaspaceandtimediscretevariationalprinciple.

By taking variations of thisspace(–time) discrete variational principle, a system of nonlinear ordinary-differential (or algebraic equations) for the free surface height and the velocity potential is obtained. Besides directly discretizing the variationalprinciple,thecentralkeyfeaturesthatmakethisdiscretizationrobustandstableareasfollows:

(i) Variationsofthepositionsofinteriornodesofthedynamicmeshmustberelatedvariationallytothefreesurfacenodes. Ifoneremovesthisintrinsicdependencetheresultingnumericalschemeispronetobeunstable.

(ii) Thenewvariationaltimediscretizationdealswiththeintrinsicallycoupleddynamicalsystem,bothintheautonomous andnon-autonomouscases,thelatterwithdrivenwaves.

Thesetwofeaturesappeartocontrastwithearlier(DG)ﬁniteelementapproaches,andresultinastablenumericalscheme withnodriftofthediscreteenergyandpreservationofdiscretephasespace.

Fora certainlinear approximationoftheflow fieldsintime, ourspace–timevariationalfiniteelement methodresults inatimediscretizationsimilartothesymplecticStörmer–Verlettimediscretization.Intheautonomouscasewithoutwave maker, we show that it is a Störmer–Verlet time discretization [15] of the dynamics atthe free surface. The latter dy-namics constituteastandardyetnonlinearcanonicalHamiltoniansystem(ofordinarydifferentialequationsbeforetimeis discretizedaswell).Thefirsttwostepsinthetimeintegrationschemeareimplicit,andthethirdstepisexplicit.Hence,the

(3)

Fig. 1. Sketch of the numerical wave basin with water waves generated by a piston wave maker.

variationalﬁniteelementdiscretizationleadstoasystemofnonlinearalgebraicequations.ItissolvedbyaNewtonmethod yieldingasparsematrix.Thesparsitydependsdirectlyonthenumberofnodessurroundingeachnodeintheﬂowdomain. Thesparsematrixstorageroutinesandsparseiterativeordirectsolversofthe PETSc packageareused

[34–36]

.

Ourvariationalﬁnite elementmethodis veriﬁedagainst anonlinear semi-analyticalsolutionby FentonandRienecker

[8]forhorizontallyperiodicwavesandvalidatedagainstexperimentaldatasetsprovidedbytheMaritimeResearchInstitute Netherlands (MARIN).Theseexperimentsare conductedusingtwo wave basin conﬁgurations,one withavarying bottom andone withaﬂatbottom.Inbothcases,wavesaregeneratedwithapistontypewave makeratoneendofthedomain. Thisexperimentalsetupcanleadtohighlyirregularwavesaswellaswavefocussing,inwhichaslowmovingwavetrainis overtakenbyafastmovingwavetrain.Moreinformationontheexperimentscanbefoundin

[16]

andanothernumerical validationinLatifahandVanGroesen

[23]

.

The outline ofthe paperis asfollows. InSection

2

, we presentMiles’ variationalprinciple for nonlinear free surface wavesandderivethegoverningequations.InSection

3

,thevariationalﬁniteelementmethodisformulated.Subsequently, the numericalveriﬁcation andexperimental validation are discussed in Section 4.Finally, conclusions are drawn in Sec-tion5.

2. Variationaldescriptionofnonlinearinviscidwaterwaves

Weconsiderwaterwavesgeneratedbyawavemakerinatwo-dimensionalwavebasininaverticalplane.Thisrestriction to two dimensionsis notessential to our approach.The ﬂuid is assumedto be inviscid,irrotational and incompressible. A sketchofthewavebasinisshownin

Fig. 1

.

Thedynamicsofwaterwavesaregovernedbythefollowingvariationalprinciple(cf.,Miles

[30]

):

0

= δ

L

(φ,

η

, φ

s

)

dt

:= δ

T

0

L R(t)

φ

s

∂

η

∂

t

−

1 2g

η

2

−

H2

dx

−

L

R(t) η

(x,t) −H(x) 1 2

|∇φ|

2_dzdx

₋

η(

R(t),t) −H(R(t)) dR dt

φ

wdz

dt

,

(1)

where

(

x

,

z

)

are the horizontal and vertical coordinates, respectively,

∇

is the gradient operator in the vertical plane,

φ (

x

,

z

,

t

)

the velocity potential, u

:= ∇φ

the irrotational velocity ﬁeld, and

Ω

= {

R

(

t

)

<

x

<

L

,

−

H

(

x

)

<

z

<

η

(

x

,

t

)

}

⊂ R

2

thetime dependentﬂowdomain.Theboundariesof

Ω

are thefreesurface boundary SF

:

z

=

η

(

x

,

t

)

with

η

thefree sur-face height, the ﬁxed bed at the bottom SB

:

z

= −

H

(

x

)

(z

=

0 is the ﬁxed referencelevel ofthe free surface at rest), a piston wave maker SW

:

x

=

R

(

t

)

and a solid wall SL

:

x

=

L. For brevity’s sake, we deﬁne

φ

s

(

x

,

t

)

= φ(

x

,

η

(

x

,

t

),

t

)

,

φ

w

(

z

,

t

)

= φ(

R

(

t

),

z

,

t

)

,

η

w

(

t

)

=

η

(

R

(

t

),

t

)

, Hw

=

H

(

R

(

t

))

and

φ

sw

(

t

)

= φ(

R

(

t

),

η

(

R

(

t

),

t

),

t

)

. Note that the wave maker movesonlyhorizontallyandthatthefreesurfaceisassumedtobesingly-valued.

Takingthevariationalderivativeswithrespectto

φ

,

η

and

φ

s in

(1)

yields

0

= δ

L

=

T

0

−

Ω(t)

∇φ · ∇δφ

dxdz

−

η

w(t) −Hw(t) dR dt

δφ

wdz

+

L

R(t)

∂

η

∂

t

δφ

s

+ φ

s

∂δ

η

∂

t

−

g

η

+

1 2

|∇φ|

2 s

δ

η

dx

−

dR dt

φ

sw

δ

η

w

dt

.

(2)

Wesimplytookthevariationswithrespecttothetime-dependentboundary z

=

η

ofthedomain

Ω(

t

)

,butcould alterna-tivelyhaveusedavariationalanalogueoftheReynoldstransporttheoremasin

[12]

.Integrating

(2)

bypartsinspaceusing Gauss’divergencetheoremaswellasintegratingtheterm

φ

s

∂δ

η

/∂

t in

(2)

bypartsintimeusingLeibnitz’s ruleyields

(4)

0

= δ

L

=

T

0

Ω(t)

∇

2

_φδφ

_dxdz

₋

∂Ω(t) n

· ∇φδφ

dS

−

η

w(t) −Hw(t) dR dt

δφ

wdz

+

L

R(t)

∂

η

∂

t

δφ

s

−

∂φ

s

∂

t

+

g

η

+

1 2

|∇φ|

2 s

δ

η

dx

+

L

R(t)

φ

s

δ

η

dx

|

0T

dt

,

(3)

with

∂Ω

:=

SF

∪

SB

∪

SW

∪

SL andn the outward unit normalvector attheboundary

∂Ω

, anddS aninﬁnitesimal line elementalongsuchaboundary.Duetotheintegrationbypartsintimeanadditionaltermarisessuchthatthelasttermin

(2),thepointatthefreesurfaceandwave-maker,cancels.Byusingtherelations

δφ

s

= (δφ)

s

+ (∂φ/∂

z

)

s

δ

η

,

(4a)

SF dS

=

L

R(t)

∇

(

z

−

η

)

dx

,

(4b) n

|

SF

= ∇(

z

−

η

)/

∇

(

z

−

η

)

,

(4c)

δ

η

(

x

,

0

)

= δ

η

(

x

,

T

)

=

0 (4d)

andaftersplittingtheboundaryintegralsoverthedifferenttypesofboundaryfaces,thevariationsin

(3)

become

0

= δ

L

=

T

0

Ω(t)

∇

2

_φδφ

_dxdz

₋

SB∪SL n

· ∇φδφ

dS

−

η

w(t) −Hw(t)

dR dt

−

∂φ

∂

x

w

δφ

wdz

+

L

R(t)

∂

η

∂

t

+

∂φ

∂

x

s

∂

η

∂

x

−

∂φ

∂

z

s

δφ

s

−

∂φ

s

∂

t

+

g

η

+

1 2

∂φ

s

∂

x

2

−

1 2

∂φ

∂

z

2 s

1

+

∂

η

∂

x

2

δ

η

dx

dt

.

(5)

Giventhearbitrarinessofthevariationsinvolvedin

(5)

,thefollowingdimensionalequationsofmotionemerge

∇

2

_φ

₌

_{0 in}

_Ω(

_t

_),

_(6a) dR dt

=

∂φ

∂

x at SW

(

t

)

:

x

=

R

(

t

),

(6b) n

· ∇φ =

0 at SB

:

z

= −

H

(

x

)

∪

SL

:

x

=

L

,

(6c)

∂

η

∂

t

+

∂φ

s

∂

x

∂

η

∂

x

−

∂φ

∂

z

s

1

+

∂

η

∂

x

2

=

0 at SF

(

t

)

:

z

=

η

(

x

,

t

),

(6d)

∂φ

s

∂

t

+

1 2

∂φ

s

∂

x

2

+

g

η

−

1 2

∂φ

∂

z

2 s

1

+

∂

η

∂

x

2

=

0 at SF

(

t

)

:

z

=

η

(

x

,

t

),

(6e)

viz. theLaplaceequation withNeumannboundaryconditions(atsolid walls)andthefree surfaceequationsexpressed in termsofthesurfacepotential

φ

s andwaveheight

η

.Recallthat

∂(φ

s

)/∂

t

= (∂φ/∂

t

)

s,cf.

(4a)

.

3. Variationalﬁniteelementformulation 3.1. Tessellation

Apoint x

∈ R

3 _in_the_space–time_domain

_E

_has_coordinates_x

_{= (}

_t

_,

_x

_,

_z

₎

_._The _ﬂow_domain

_Ω(

_t

₎

_at_time _{t is}_deﬁned_as

Ω(

t

)

:= {(

x

,

z

)

∈ R

2

|

(

t

,

x

,

z

)

∈

E}

.Thespace–timedomainboundary

∂E

consistsofthehypersurfaces

Ω(

0

)

:= {

x

∈ ∂

E |

t

=

0

}

,

Ω(

T

)

:= {

x

∈ ∂

E |

t

=

T

}

and

Q

:= {

x

∈ ∂

E |

0

<

t

<

T

}

=:

SF

∪

SW

∪

SB

∪

SL

× (

0

,

T

)

.

Usingthepartition0

≤ · · · <

tn

<

· · · ≤

tN+1

=

T ofthetimeinterval

(

0

,

T

)

,wecansubdividethespace–timedomain

E

intospace–timeslabs

E

n

_{:= {}

_x

_∈

_{E |}

_t

_{∈ (}

_tn

_,

_tn+1

₎

_}

_with₀

_≤

_n

_≤

_N. _If_the_space–time_boundary_faces _that_are_part_of

_Q

_are

restrictedtothespace–timeslab

E

nthenweusethesuperscriptn.Theremainingboundariesofthespace–timeslab

E

n are

Ω

n

= Ω(

tn

)

and

Ω

n+1

= Ω(

tn+1

)

.

The tessellation

T

h isconstructed bygenerating amesh withquadrilateralelements inthedomain

Ω(

t

)

,yielding the computational domain

Ω

h

(

t

)

.Thespatial elementsaredenoted asKke,withke theelementindexinthetessellation.Due

(5)

tothefreesurfaceandwavemakermotionthenodesofeachelement K_kn e att

=

t

n _in_domain

_Ω

h

(

tn

)

willmovetoanew positionatt

=

tn+1,resultinginthespatialelement Kn_k+1

e anddomain

Ω

h

(

t

n+1

₎

_._The_mesh_deformation_due_to_the_domain

boundarymotionisunknownaprioriandpartofthesolution.Thegeometryoftheelements isthereforeupdated during thesolutionprocess.

Wewillconsiderstructuredmeshes.Nodesatthefreesurfaceandwavemakerbelong,respectively,tothesets

N

F and

N

W; nodesinthedomain

Ω

h are collected intheset

N

Ωh andnodesofthe domainexcludingthefree surface nodesin

¯

Ω

h

\

SF,with

Ω

¯

h theclosureof

Ω

h, intheset

N

Ω¯h. Thedomain issplit intwo parts:in

Ω1h

:

R

(

t

)

<

x

<

Lw the nodes inthemeshcanmove horizontallyaswellasverticallyandin

Ω2h

:

Lw

<

x

<

L with Lw

<

L and

−

H

(

x

)

<

z

<

η

(

x

,

t

)

the nodesinthemeshmoveonlyinthevertical.

3.2. Functionspacesandapproximations

Thespacefortwo-dimensionalbasisfunctions

ϕ

˜

on

Ω

forthepotential

φ

isW1

(Ω)

= { ˜

ϕ

∈

H1

(Ω)

}

withthestandard Sobolevspace H1

_(Ω)

_._The _{approximation}

_φ

h to potential

φ

is

φ

h

∈

W kp h

=

W1

(Ω)

∪

X kp h with X kp

h a ﬁniteelement space of continuous Lagrange basis functions, and includingat least polynomials of degree kp on each element Kke. The free surface SF is parameterizedwithcoordinatex

∈ [

R

(

t

),

L

]

.The associatedspacefortheone-dimensional basisfunctions

ϕ

onx

∈ [

R

(

t

),

L

]

forthefreesurfaceheight

η

andthefreesurfacepotential

φ

s isW1

(

[

R

(

t

),

L

])

= {

ϕ

∈

H1

(

[

R

(

t

),

L

])}

.

The nodes of the elements are time dependent, making the standard global basis functions also time dependent, as follows.Forsimplicity,we usequadrilateralelements withthestandardtwo-dimensionalpiecewise linearandcontinuous (global)basis functions

ϕ

˜

j

(

x

,

z

,

t

)

inspaceforthe“interior”potential

φ (

x

,

z

,

t

)

.Thefourlocalbasisfunctionsona quadri-lateralelement Kke are

ˆψ

α

(ζ

1

, ζ

2

)

=

1

4

(

1

± ζ

1

)(

1

± ζ

2

)

(7)

on the reference element

K

= (ζ

1,

ζ2)

T

= [−

1

,

1

]

2 with local node numbers

α

=

1

,

. . . ,

4.Physical coordinates andlocal referencecoordinatesarerelatedbyamappingFke withthesamestandardshapebasisfunctions

Fke

:

K

→

Kke

: (ζ

1

, ζ

2

)

→ (

x

,

z

)

:= (

xke(α)

,

zke(α)

) ˆ

ψ

α

(ζ

1

, ζ

2

),

(8) withthe fourvertex coordinates

(

xke(α)

(

t

),

zke(α)

(

t

))

of Kke attime t, globalnode number j

=

ke

(

α

)

, inwhich Einstein’s convention ofsummation over repeatedindices(here

α

) has beenused. Similarly,one-dimensional piecewise linear and continuousbasis functions

ϕ

l

(

x

,

t

)

inspaceareusedatthefreesurface,for

η

(

x

,

t

)

and

φ

s

(

x

,

t

)

.Givenfree-surfaceelement

¯

Kks andlocalnodenumber

β

,theglobalfree-surfacenodenumberl

=

ks

(β)

.Themappingforeachfreesurfaceface K

¯

ks to thereferenceelement K

¯

= ζ ∈ [−

1

,

1

]

reads

¯

Fs

: ¯

K

→ ¯

Kks

: ζ →

x

=

xks(β)

(

t

)ψ

β

(ζ ),

(9)

withx-coordinatesxks(β)

(

t

)

ofthetwosurfacenodes.Thestandardtwoelementalshapefunctionswithks theglobal free-surface element indexare

ϕ

ks(β)

(

x

(ζ,

t

),

t

)

= ψ

β

(ζ )

=

1₂

(

1

± ζ)

withlocalindex

β

=

1

,

2. Whenreferencecoordinates are

usedthetimedependenceintheformulationappearsintheJacobianofthemapping(seealso

Appendix A

). Theresultingglobalﬁniteelementapproximationsare

φ

h

(

x

,

z

,

t

)

= φ

j

(

t

)

ϕ

˜

j

(

x

,

z

,

t

),

(10a)

φ

sh

(

x

,

t

)

= φ

k

(

t

)

ϕ

k

(

x

,

t

)

and

η

h

(

x

,

t

)

=

η

l

(

t

)

ϕ

l

(

x

,

t

),

(10b)

inwhichindicesi

,

j

=

1

,

. . . ,

Nn concernallnodesofthemeshandindicesk

,

l

,

r

,

s

=

1

,

. . . ,

Nf onlythefreesurface nodes (indicesr

,

s areusedlater).Summationoverrepeatedindicesisadopted.Hereafter,nodesindexedbyi

,

jconcernallnodes exceptthefree surfaceones. Duetothe deﬁnitioninreferencecoordinates,it followsthat

(

ϕ

˜

l

)

s

=

ϕ

l,i.e., alignedhereby taking

ζ2

=

1 (see

Appendices

AandB).

Followingourrestrictiontoasingle-valuedfreesurfaceandaparameterizationofthefreesurfaceasSF

: (x

,

z

=

η

(

x

,

t

))

, the x-coordinate of the free surface nodes is a prescribed time-dependent function xk

(

t

)

in the x-direction, while the z-coordinate is z

=

η

h

(

xk

,

t

)

. When there is no wave maker motion, the horizontal node positions xi are ﬁxed in time. A consequence of thisrestriction is that the movementsof thevertically-aligned interiornodes are slavedto the corre-spondingfreesurfacenodesattherespectivex-location.Thecomputationalmeshinthecomputationalﬂowdomain

Ω1h

(

t

)

isthenobtainedbymovingeachnodeikaccordingto:

xik

(

t

)

=

xik

(

0

)

+

γ

x,ikR

(

t

)

and zik

(

t

)

=

zik

(

0

)

+

γ

z,ik

η

(

xk

,

t

),

(11) withthe multiplication factorsgiven by

γ

x,ik

:= (

Lw

−

xik

(

0

))/

Lw and

γ

z,ik

:= (

H

(

xk

)

+

zik

(

0

))/

H

(

xk

)

. Nodes ik belowthe respective free surface node k lie on vertical mesh lines. For the free surface node k, we ﬁnd indeed that zk

=

η

(

xk

,

t

)

since zk

(

0

)

=

0. Forthe one at the bottom zik

(

t

)

= −

H

(

xk

)

since then zik

(

0

)

= −

H

(

xk

)

. When the initial free surface at t

=

0 hasno elevation,then zk

(

0

)

=

0 for nodesatthefree surface.So herewe havelimitedthenodes xik

(

0

)

and zik

(

0

)

to form a structured mesh. In

Ω

2h, only the vertical movement in (11)is retained withno horizontal node movement.

(6)

Fig. 2. AsketchismadeofthediscretewavebasinwithameshintheareaR(t)<x<Lwand

−

H(x)<z<η(x,t)inﬂuencedbythepistonwavemakerat

x

=

R(t)withbothhorizontalandverticalmeshmovement,andanareaLw<x<L withonlyverticalmeshmovement.

A corresponding sketchofthediscretizeddomainisgivenin

Fig. 2

.Ourﬁniteelementapproachpermitsthislimitation to beliftedforageneralparameterizationSF

: (

xs

(

t

),

zs

(

t

))

givenfreesurfacenodesk withpositions

(

xs,k

,

zs,k

)

.

Theconsequenceofthemeshmovementisthereforethatthebasisfunctionsin

(10)

dependexplicitlyontime.Evenfor thecaseofunstructuredmeshmovement,thistimedependencecangenericallybeexpressedas

˜

ϕ

j

(

x

,

z

,

t

)

= ˜

ϕ

j

x

,

z

;

η

,

R

(

t

)

and

ϕ

j

(

x

,

t

)

=

ϕ

j

x

;

R

(

t

)

.

(12)

The vector

η

(

t

)

has components

η

r

(

t

)

, the vector

˜φ(t

)

components

φ

i

(

t

)

at interior nodes, and

φ

components

φ

l

(

t

)

at free surface nodes.Hence, basisfunctionsimplicitlydependontime throughthetime dependenceofthe variables

η

and explicitly depend on time through the prescribed function R

(

t

)

for the piston position. The basis functionson the free surfacedothereforenotdependontimewhenR

(

t

)

isconstant,i.e.,intheabsenceofawavemaker.Hereafterwemaynot alwaysdistinguishthesedifferenttypesoftimedependence inthebasisfunctionsexplicitly.

3.3. Space-plus-timeformulation

Wedemonstratehowthetwokeychallenges(i)and(ii)mentionedinSection1areresolvedbyﬁrstdiscretizingthe vari-ationalprinciple

(1)

inspaceandsubsequentlyintime.WithoutanyintegrationbypartsoruseofGauss’law,substitution ofexpansions

(10)

intothevariationalprinciple

(1)

yields

0

= δ

T

0 L

˜

φ, φ,

η

;

R

(

t

)

dt

= δ

T

0

Mkl

φ

k d

η

l dt

−

Dkl

φ

k

η

l

−

1 2g Mkl

η

k

η

l

−

1 2Ai j

φ

i

φ

j

−

Wm

φ

m

dt

,

(13)

withthefollowingmatricesandvectorsdeﬁnedas

Mkl

R

(

t

)

:=

L

R(t)

ϕ

k

ϕ

ldx

,

Dkl

R

(

t

)

:= −

L

R(t)

∂

ϕ

l

∂

t

ϕ

kdx

,

(14a) Ai j

η

(

t

),

R

(

t

)

:=

Ωh(t)

(

∇ ˜

ϕ

i

· ∇ ˜

ϕ

j

)

dxdz

,

(14b) Wm

η

1

(

t

),

R

(

t

)

:=

η

1(t) −Hw(t) dR dt

ϕ

˜

m

x=R(t) dz

,

(14c)

where

η

1(t

)

isthewaveheightvariableonthewavemaker,andindexm isusedforthenodesonthewavemaker.Theﬁrst node withm

=

1 liesatthefreesurfaceandallothernodesonthewavemakercarryindexm.Wewillalsouseindexm

˜

fortheseinteriorwavemakernodes.Notethatitmaybe desirabletodiscretizethebottomboundaryalsoinapiecewise linearmanner.Nearthewavemakerthebottomis,forsimplicity,takentobeﬂat:Hw

=

H

(

R

(

t

))

isconstant.Thesedifferent kindsoftime dependencehavebeenexplicitlyindicatedin

(14)

followingthatdependencein

(12)

.ThedependenceAi j

(

η

)

signiﬁesthat Ai j inprincipledependsonall

η

k’s.

Variationof

(13)

whileincorporating

(14)

yields

0

=

T

0

Mkl d

η

l dt

−

Dkl

η

l

−

Aik

φ

i

−

W1

δ

1k

δφ

k

−

d

(

Mkl

φ

k

)

dt

+

Dkl

φ

k

+

g Mkl

η

k

+

1 2

∂

Ai j

∂

η

l

φ

i

φ

j

+

∂

Wm

∂

η

1

δ

1l

φ

m

δ

η

l

− (

Aij

φ

i

+

Alj

φ

l

+

Wm

δ

mj

)δφ

j

dt

,

(15)

(7)

inwhichtheboundarytermarisingfromtheintegrationbypartsintimehasvanishedbecauseweused

δ

η

l

(

0

)

= δ

η

l

(

T

)

=

0 (the ﬁrst relation says that the variation of the initial condition is zero, while the latter follows similarly by reversing time),andinwhich

δ

mk istheKroneckerdeltasymbol:itisonewhenm

=

k andzerootherwise.Theordinarydifferential equationsaretherefore

δ

η

l

:

d

(

Mkl

φ

k

)

dt

+

Dkl

φ

k

+

g Mkl

η

k

+

1 2

∂

Ai j

∂

η

l

φ

i

φ

j

+

∂

Wm

∂

η

1

δ

1l

φ

m

=

0 (16a)

δφ

k

:

Mkl d

η

l dt

−

Dkl

η

l

−

Aik

φ

i

−

W1

δ

1k

=

0 (16b)

δφ

j

:

Aij

φ

i

+

Alj

φ

l

+

Wm

δ

mj

=

0

,

(16c)

withclearlydistinguishedvariationsofsurfacedegreesoffreedomarisingfromthearbitrarinessof

δ

η

l and

δφ

k,andinterior onesarisingfrom

δφ

j.

ThediscretizedLaplaceequationrevealsthatwecaneliminatetheinteriordegreesoffreedomindexedbyi

,

jasfollows

φ

i

= −

AljA−_j_i1

φ

l

−

WmA−_m1_i

.

(17) Substitution of (17)into the variationalprinciple

(13)

then leads to Hamiltonianfree surface dynamicsgoverned by the followingvariationalprinciple

0

= δ

T

0 L

φ,

η

;

R

(

t

)

dt

= δ

T

0

Mkl

φ

k d

η

l dt

−

1 2g Mkl

η

k

η

l

−

1 2Bkl

φ

k

φ

l

−

Dkl

φ

k

η

l

−

Cl

φ

l

−

Fw

dt

,

(18a)

inwhichtheSchurcomplement

Bkl

η

,

R

(

t

)

=

Akl

−

AkiA−_i1jAjl

(18b)

ofthefull matrix Ai j emerges asa consequenceofthiseliminationoftheinteriordegreesoffreedom,a vector Cl anda function Fw,i.e., Cl

η

,

R

(

t

)

=

W1

δ

1l

−

AljA−_j1_mWm and Fw

(

t

)

= −

1 2Wm˜A −1 ˜ mmWm

.

(18c)

Thenotation Bkl

(

η

)

meansthat Bkldependsinprincipleonallsurfacedegreesoffreedom

η

.

Variationalprinciple

(18)

inessencecapturesadiscreteboundaryelementrepresentationofthedynamics.Onecanverify thatthedynamics

(16)

,inwhichtheinteriordegreesoffreedom

φ

i or

φ

j havebeeneliminatedusing

(17)

,followsdirectly fromthisfree-surfacevariational principle

(18)

.Wenevercalculatethematrixinverse A−_i1j explicitly.Thisinverseisonly usedtoderivetheappropriatetimesteppingschemeinthenextstep.Intheend,we“unfold”thealgebraicsystemagainto avoidadirectcalculationofthisinversematrixbecausethatiscomputationallyexpensive.Itislargeanddependsimplicitly andexplicitlyontime.

Thevariationalprinciple

(18)

isacanonicalHamiltoniansystem

0

= δ

T

0

pl d

η

l dt

−

H

p,

η

;

R

(

t

)

dt (19)

withconjugatevariables

{

pl

,

η

l

}

= {

pl

=

Mkl

φ

k

,

η

l

}

andextended(Hamiltonian)functionH

=

Hh

+

Hna splitintothe follow-ingHamiltonian, Hh,andwavemakerpart,Hna (non-autonomous),

Hh

≡

1

2g Mkl

η

k

η

l

+

1

2Bkl

φ

k

φ

l and Hna

≡

Dkl

φ

k

η

l

+

Cl

φ

l

+

Fw

.

(20)

Without wave maker Hh is time independent and Hna disappears. A standard yet semi-implicit and symplectic ﬁnite-differenceStörmer–Verlettimeintegrator

[15]

wouldthensuﬃcetodiscretize

(18)

forR

(

t

)

=

0 (or

(19)

withHna

=

0).

The explicit time dependence that makes (18) a non-autonomous variational principle stimulated us to develop a discontinuous Galerkin ﬁnite element approach in time. In each time slab the DG basis functions

ψ

t

(

t

)

are continuous and the associated space is Vkt

h

= {ψ

t

∈

P

kt

(

t

∈ (

t

n

_,

_tn+1

₎₎

_}

_with

_P

kt polynomials of degree kt. We will restrict atten-tion to piecewise linearpolynomials

ψ

t

(

t

)

,i.e., kt

=

1. Wechoose the followinglinear expansionswithin each time slab t

∈ (

tn

_,

_tn+1

₌

_tn

₊

_tn

₎

_:

(8)

Fig. 3. Piecewise linear approximation in time of the potentialφ, free surface potentialφsand wave heightη.

η

l

(

t

)

=

η

n_l,+

tn+1

−

t

tn

+

η

n_l+1,−

t

−

tn

and (21a)

φ

k

(

t

)

= φ

_kn+1/2 2

(

t

−

tn

)

tn

+ φ

n,+ k

(

tn

+

tn+1

−

2t

)

tn

,

(21b)

where

η

_ln,±

=

lim→0

η

l

(

tn

±

)

,

φ

kn,+

=

lim→0

φ

k

(

tn

+

)

and

φ

n+1/2

k

= φ

k

(

t

=

tn+1/2

)

.(For symplectic schemes thetime step needs to be ﬁxed, such that

tn

=

t, also in the simulations presented later.) A sketch of these temporal basis functionsisgivenin

Fig. 3

.Theseexpansionsarecontinuouswithinthetimeslabbutdiscontinuousacross.Hence,wehave todeﬁnethedeltafunctionarisingfromthetimederivativeinthetermMkl

φ

kd

η

l

/

dt in

(18)

.Thecontributionofthisterm can eitherbeposited, introducedbyintegratingbypartstwicewhileintroducinganumericalﬂux inonestepbutnotthe other,orbyappealingtothetheoryofnonconservativeproducts(forhyperbolicsystems

[7,32,31]

orquasi-linearequations

[41]).Alternatively,forthevariationalprinciplesconsideredthiscontributionarisesasalimitofacontinuousGalerkinﬁnite elementmethodintimewithevenandoddnumberedelementsinwhichtheoddelementsarelimitedtozero.1 _In_all_these

differenttypesofcalculations,atnodetn+1 anoutcomeisthecontribution:

Mkl

φ

k d

η

l dt

≈

M n+1,+ kl

φ

n+1,+ k

η

_ln+1,+

−

η

_ln+1,−

,

(22)

i.e., Mkl

φ

k is evaluated on the future side of node tn+1 anda jumpin

η

across tn+1 emerges. When Mkl is considered continuousintime,i.e.,forcontinuous R

(

t

)

,then M_kln+1,+

=

Mn_kl+1.Wewilltake R

(

t

)

continuous.

The other termsinthetime integralof

(18)

are approximatedovereach time slabsuch thatthey aresecond orderin time.Theresultingintegralsarethenapproximatedbythefollowingquadraturerules:

tn+1

tn Mkl

φ

k d

η

l dt dt

≈

M n+1/2 kl

φ

n+1/2 k

η

n_l+1,−

−

η

_ln,+

;

(23a) tn+1

tn Dkl

φ

k

η

ldt

≈

1 2

t n_Dn+1/2 kl

φ

n+1/2 k

η

n_l,+

+

η

n_l+1,−

;

(23b) tn+1

tn 1 2Mkl

η

k

η

ldt

≈

1 4

t n_Mn+1/2 kl

η

_kn,+

η

n_l,+

+

η

_kn+1,−

η

n_l+1,−

;

(23c) and, tn+1

tn 1 2Bkl

φ

k

φ

ldt

≈

1 4

t n

_Bn+1/2 kl

η

n,+

₊

_Bn+1/2 kl

η

n+1,−

_φ

n+1/2 k

φ

n+1/2 l

,

(23d) 1 _See_{Appendix C.}

(9)

tn+1

tn Cl

φ

ldt

≈

1 2

t n

_Cn+1/2 l

η

n,+

₊

_Cn+1/2 l

η

n+1,−

_φ

n+1/2 l

.

(23e)

Intheﬁrstintegral,weeitherapproximateMkl atthemid-timelevel,orwecanviewitasanexactevaluationoftheentire conjugatevariableexpandedinalinearwayintime.Theremainingintegralsareapproximatedusingamid-pointrulewith

φ

k and time evaluated at t

=

tn+1/2.Here, C n+1/2

l

(

η

n,+

)

=

Cl

(

η

n,+

,

R

(

tn+1/2

))

, etc. The time derivative of

ϕ

in Dkl leads via thechainrule toa partialderivative withrespectto R

(

t

)

timesdR

/

dt.The latteriseither foundexplicitlyor R

(

t

)

is approximatedandthenallexpressions areevaluated attn+1/2.Hence,bycombiningtheresultingapproximationsderived in

(22)

and

(23)

thealgebraic,space–timediscretevariationalprincipleatthefreesurfaceinessencebecomes

0

= δ

N

n=0 L

φ

n,+

, φ

n+1/2

,

η

n,+

_,

_η

n+1,−

_(24a)

= δ

N

n=0 M_kln+1/2

φ

_kn+1/2

η

_ln+1,−

−

η

n_l,+

−

t n 2

Hn+1/2

φ

n+1/2

,

η

n,+

₊

_Hn+1/2

_φ

n+1/2

_,

_η

n+1,−

+ δ

N

n=−1 M_kln+1

φ

_kn+1,+

η

_ln+1,+

−

η

n_l+1,−

,

(24b) withHn+1/2

_(φ

n+1/2

,

η

n+1,−

₎

₌

_H

_(φ

n+1/2

,

η

n+1,−

_,

_tn+1/2

₎

_.

Variation of(24) withrespect to the independentvariables indicated yields an extension of the symplectic Störmer– Verlettimeintegrationscheme

δφ

_kn+1,+

:

η

_ln+1,+

=

η

n_l+1,− (25a)

δ

η

_ln,+

:

M_kln+1/2

φ

_kn+1/2

=

Mn_kl

φ

_kn,+

−

1 2

t n

∂

Hn+1/2

(φ

n+1/2

,

η

n,+

)

∂

η

_ln,+ (25b)

δφ

_kn+1/2

:

M_kln+1/2

η

_ln+1,−

=

M_kln+1/2

η

n_l,+

+

1 2

t n

∂

Hn+1/2

(φ

n+1/2

,

η

n,+

₎

∂φ

n_k+1/2

+

∂

Hn+1/2

(φ

n+1/2

,

η

n+1,−

₎

∂φ

_kn+1/2

(25c)

δ

η

_ln+1,−

:

M_kln+1

φ

_kn+1,+

=

M_kln+1/2

φ

_kn+1/2

−

1 2

t n

∂

Hn+1/2

(φ

n+1/2

,

η

n+1,−

)

∂

η

_ln+1,−

.

(25d)

Weobservethatwhiletheﬁniteelementexpansionintimewasoriginallyposedtobediscontinuous,thevariationsenforce

η

l to be continuous while

φ

k remains discontinuous in time. Both the choice of the numerical ﬂux (or path) and the quadraturesusedabovewereguidedbyourintentiontoderivetheclassicStörmer–Verletschemeintheabsenceofexplicit timedependence,i.e.,herewhenR

(

t

)

=

0.

Finally,theinverseofthematrix Aij inthenon-constantmatrixBklandthenon-constantvector Cl canbeavoidedby unfoldingthefreesurfacedynamicsagaintoincludetheinteriordegreesoffreedom.Whenweuse

φ

_in,+

= −

A_ljn+1/2

η

n,+

A−_j_i1

η

n,+

,

tn+1/2

φ

_ln+1/2 (26a)

φ

n_i+1,−

= −

Anlj+1/2

η

n+1,−

_A−1 ji

η

n+1,−

_,

_tn+1/2

_φ

n+1/2 l (26b)

attwotimelevels,theintegralofthekineticenergy K inthestandardHamiltonianHh

=

K

+

P

=

Ai j

φ

i

φ

j

/

2

+

g Mkl

η

k

η

l

/

2 canberewrittenas2 N

n=0 1 2

t n

_K

˜

_φ

n,+

_{, φ}

n+1/2

_,

_η

n,+

_,

_tn+1/2

₊

_K

˜

_φ

n+1,−

_{, φ}

n+1/2

_,

_η

n+1,−

_,

_tn+1/2

≡

N

n=0

tn 4

Bn_kl+1/2

η

n,+

₊

_Bn+1/2 kl

η

n+1,−

_φ

n+1/2 k

φ

n+1/2 l

=

N

n=0

tn 4

A_kln+1/2

η

n,+

φ

_kn+1/2

φ

n_l+1/2

+

An_kl+1/2

η

n+1,−

φ

n_k+1/2

φ

_ln+1/2

2 _In_(18a)_we_derived_that_A

i jφiφj=Bklφkφlusing(17)with(18b).Herein(27),startingfrom(23d),wetracethestepsinthisderivationbackbyusing

(10)

+

An_kj+1/2

η

n,+

_φ

n+1/2 k

φ

n,+ j

+

A n+1/2 kj

η

n+1,−

_φ

n+1/2 k

φ

n+1,− j

+

An_i+l1/2

η

n,+

_φ

n,+ i

φ

n+1/2 l

+

A n+1/2 il

η

n+1,−

_φ

n+1,− i

φ

n+1/2 l

+

An_i+_j1/2

η

n,+

φ

_in,+

φ

n_j,+

+

An_i+_j1/2

η

n+1,−

φ

_in+1,−

φ

n_j+1,−

(27) after some algebraic manipulations. A similar unfolding can be done for the wave maker term involving Cl but is not required.Hence,itseemsthattheﬁniteelementexpansionfortheinteriordegreesoffreedomis

φ

i

= φ

n_i,+

tn+1

−

t

tn

+ φ

n+1,− i

t

−

tn

(28) with

φ

_in,±

=

lim→0

φ

i

(

tn

±

)

.

Analogously,wenowseethatthetimeintegralofthewavemakercontributionin

(13)

isapproximatedas T

0 Wm

φ

mdt

∼

=

N

n=0

tn 2

Wqn+1/2

η

n,+ q

φ

qn+1/2

+

Wqn+1/2

η

n+1,− q

φ

qn+1/2

+

W_mn+1/2

η

n,+ q

φ

n_m,+

+

Wmn+1/2

η

n+1,− q

φ

n_m+1,−

,

(29)

withinteriorwavemakernodesm,thewavemakernodeq

=

1 atthefreesurfaceandW_mn+1/2

(

η

nq,+

)

=

Wm

(

η

nq,+

,

tn+1/2

)

, etc.Theverticesonthewavemakerarethusevaluatedattimetn+1/2.

Byusingapproximations

(23a)–(23c)

fortheﬁrstthreefreesurface andpotentialenergytermsin

(13)

and approxima-tions

(27)

forthekineticenergyand

(29)

forthewavemakertermsin

(13)

,thespacediscreteprinciple

(13)

becomes the followingspace–timediscretevariationalprinciple

0

= δ

L

˜

φ

n,+

, ˜

φ

n−1,−

, φ

n,+

, φ

n+1/2

,

η

n,+

,

η

n+1,−

= δ

N

n=0 Mn_kl+1/2

φ

_ln+1/2

η

_kn+1,−

−

η

n_k,+

−

1 2

t n_Dn+1/2 kl

φ

n+1/2 l

η

_kn,+

+

η

_kn+1,−

−

1 2

t n

_Hn+1/2 h

˜

φ

n,+

, φ

n+1/2

,

η

n,+

₊

_Hn+1/2 h

˜

φ

n+1,−

, φ

n+1/2

,

η

n+1,−

−

1 2

t n

_Wn+1/2 q

η

n,+ q

+

Wqn+1/2

η

n+1,− q

φ

qn+1/2

+

Wmn+1/2

η

n,+ q

φ

_mn+,

+

Wmn+1/2

η

n+1,− q

φ

_mn+1,−

+

N

n=−1 M_kln+1

φ

_ln+1,+

η

_kn+1,+

−

η

_kn+1,−

,

(30)

withindicesk

,

l

∈

N

F,m

∈

N

W

\

N

F,i

,

j

∈

N

Ω¯h andq

∈

N

F

∩

N

W.Itconsistsofsumsofdiscretizedintegralsandjump termsoverandacrossthespace–timeslabs

E

n.

3.4. Dynamics

Takingvariationsof

(30)

overtheindependentvariables

φ

_ln+1/2

,

φ

_ln,+,

φ

n_i,+,

φ

in+1,−,

η

n,+

k and

η

n+1,−

k yieldsanalgebraic systemforthesepotentialandwaveheightcoeﬃcients:

δφ

_ln+1,+

:

M_kln+1

η

_kn+1,+

=

Mn_kl+1

η

n_k+1,−

,

(31a)

δφ

_in,+

:

∂

Hn_h+1/2

( ˜

φ

n,+

, φ

n+1/2

,

η

n,+

₎

∂φ

n_i,+

= −

W_mn+1/2

η

nq,+

δ

mi

,

(31b)

δ

η

n_k,+

:

Mn_kl+1/2

+

1 2

t n_Dn+1/2 kl

φ

_ln+1/2

=

Mn_kl

φ

n_l,+

−

1 2

t n

_∂

_Hn+1/2 h

( ˜

φ

n,+

, φ

n+1/2

,

η

n,+

₎

∂

η

n_k,+

+

∂

Wnq+1/2

(

η

nq,+

)

∂

η

nq,+

φ

qn+1/2

δ

kq

+

∂

W_mn+1/2

(

η

qn,+

)

∂

η

qn,+

φ

_mn,+

δ

kq

,

(31c)

δφ

_ln+1/2

:

Mn_kl+1/2

−

1 2

t n_Dn+1/2 kl

η

n_k+1,−

=

Mn_kl+1/2

+

1 2

t n_Dn+1/2 kl

η

_kn,+

+

1 2

t n

_∂

_Hn+1/2 h

( ˜

φ

n,+

, φ

n+1/2

,

η

n,+

₎

∂φ

_ln+1/2

+

∂

Hn_h+1/2

( ˜

φ

n+1,−

, φ

n+1/2

,

η

n+1,−

₎

∂φ

_ln+1/2

Variational space–time (dis)continuous Galerkin method for nonlinear free surface water waves

Journal

of

Computational

Physics

Variational

space–time

(dis)continuous

Galerkin

method

for

nonlinear

free

surface

water

waves

E. Gagarina

,

∗

,

V.R. Ambati

,

J.J.W. van der Vegt

,

O. Bokhove

,

,

∗

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

[18]

*

[26]

[30]

[21]

[25]

[40]

[33]

[37]

[3,5,9,10,14,17]

[2]

+

[38]

[34–36]

[16]

[23]

2

3

Fig. 1

[30]

= δ

L

(φ,

η

, φ

)

:= δ

φ

∂

η

∂

−



₋