Time-dependent conformal mapping of doubly-connected regions

This paper examines two key features of time-dependent conformal mappings in doubly-connected regions, the evolution of the conformal modulus Q(t) and the boundary transformation generalizing the Hilbert transform. It also applies the theory to an unsteady free surface flow. Focusing on inviscid, incompressible, irrotational fluid sloshing in a rectangular vessel, it is shown that the explicit calculation of the conformal modulus is essential to correctly predict features of the flow. Results are also presented for fully dynamic simulations which use a time-dependent conformal mapping and the Garrick generalization of the Hilbert transform to map the physical domain to a time-dependent rectangle in the computational domain. The results of this new approach are compared to the complementary numerical scheme of Frandsen (J. Comput. Phys. 196:53–87, 14) and it is shown that correct calculation of the conformal modulus is essential in order to obtain agreement between the two methods.


Introduction
Conformal mapping of a simply connected region, with a suitably smooth boundary, to a fixed disc, in the steady or time-dependent case, is assured by the Riemann mapping theorem. The Riemann map becomes unique when three points of the map are fixed. The interior of the region can be either mapped to the interior of the disc or to the periodic half space as shown in Fig. 1. Since the application of interest in this paper is water wave problems it is the mapping to the periodic half space that is of interest. One of the great advantages from a numerical point of view is that the real and imaginary parts of a complex analytic function restricted to the boundary can be related by the Hilbert transform, which is amenable to fast numerical evaluation.
When the region is doubly connected, it can be mapped to an annulus or to the infinite periodic strip, as shown in Fig. 2. Two new significant problems arise in the doubly connected case. The first problem is that doubly-connected regions have a conformal invariant called the conformal modulus. The mapping is conformal only for certain values of the conformal modulus (e.g. Chapter 17 of [19,25]). If the mapping is time-dependent then in general the conformal modulus will also be time dependent.
The second problem that arises is that the boundary values of the real and imaginary parts of a complex analytic function are no longer simply related by the Hilbert transform. Indeed, the boundary values on the inner radius are intrinsically connected to the boundary values on the outer radius. The modified boundary transform was first worked out by Garrick [16] for the annulus, and a detailed derivation is given in Section 17.4 of Henrici [19] and Section 1.8 of Papamichael and Stylianopoulos [25]. The modification of the Hilbert-Garrick transform to the periodic strip is straightforward in principle, but we will sketch the derivation here as there are some subtle sign changes.
In this paper we study numerically the dynamics of the conformal modulus associated with a time-dependent conformal mapping, derive the Hilbert-Garrick transform, and apply them to fluid sloshing of a finite-depth fluid. A typical configuration of interest is shown in Fig. 3. The walls x = 0, x = 1 and y = 0 are fixed and the upper boundary, defined parametrically by x = X(μ, t), y = Y (μ, t), for 0 ≤ μ ≤ 1, is a time-dependent curve. The image of this domain under a conformal map is the time-dependent region  where Q(t) is the time-dependent conformal modulus. The motivation is the two-dimensional space and time periodic water wave problem in a bounded domain, otherwise known as the "sloshing" problem [13]. The shape of the fluid domain changes with time as the free-surface evolves, and the position of the free-surface is not known a priori. There is a vast array of numerical methods for time-dependent free-surface problems (cf. [10]). One approach is to transform the moving domain to a fixed domain, in particular to a rectangular or annular domain, where analytic and computational approaches can be more readily applied. Of interest here are mappings that are conformal. Time-dependent conformal mappings have been widely used in the simulation of water waves in the case of infinite depth where a conformal modulus is not needed (e.g. [4-7, 11, 17, 23, 24, 26, 27, 31, 32] and references therein). An additional simplification in infinite depth is the use of the Hilbert transform to connect the real and imaginary parts of the mapping at the boundary.
Conformal mapping is most effective for small to moderate amplitude waves, and this is the regime that we will concentrate on in this paper. For large amplitude waves and sharp interfaces conformal mapping is less effective (grid points at the surface are badly distributed near sharp crests) and methods such as the boundary element method are more effective (e.g. [33]).
In finite depth, conformal mapping for water wave dynamics brings in two difficulties: the appearance of the conformal modulus and the need to extend the Hilbert transform to the Hilbert-Garrick transform. Whereas the Hilbert transform is amenable to fast transform methods (e.g. [18]), the Hilbert-Garrick transform depends explicitly on the conformal modulus and so has to be recalculated at every Fig. 3 Schematic of the conformal mapping of the time-dependent region to a time-dependent rectangle update of the conformal modulus. Finite depth simulations using conformal mappings have been reported in a number of papers (e.g. [6,8,12,22,30]). However, there does not appear to be a systematic study of the implications of a time-dependent conformal modulus.
Unfortunately, the conformal modulus gives an additional unknown that has to be computed. It is a purely geometric requirement of the mapping and does not appear to have any physical significance. Hence there is a temptation to find a strategy that allows fixing the conformal modulus at a chosen value. When the flow is steady the mapping is fixed and the conformal modulus is fixed. In this case, if the specific value of the mean depth is not important, the value of the conformal modulus can be fixed, and then the value of the depth determined as part of the solution. This strategy is used in Vanden-Broeck and Schwartz [29] and Constantin et al. [9]. However, this strategy fails in the time-dependent case since fixing the conformal modulus will lead to a time-dependent value of the depth! In this paper time-dependent conformal mappings are considered such that the physical Cartesian coordinates (x, y) transform as where (μ, ν) are coordinates for V (t) in Eq. 1.1 and t is time. These coordinates satisfy the Cauchy-Riemann equations x μ = y ν and x ν = −y μ , (1.3) in the mapped domain, and thus both satisfy Laplace's equation A kinematic approach to studying the time-dependence of the conformal modulus is given in Section 2, where a time-dependent free surface evolution is explicitly specified, as a model of linear free sloshing in a rectangular vessel. It is shown analytically by applying small time asymptotics that the conformal modulus is indeed time-dependent, and it is concluded that if it is not calculated correctly, or chosen to be a fixed constant, then O(1) errors may occur in the solution.
In Section 3 the generalization of the Hilbert transform to the Hilbert-Garrick transform is sketched for the periodic infinite strip. In Section 4, the time-dependent conformal mapping of the sloshing problem with an unknown free surface is computed with time-dependent modulus. This numerical approach is appealing over schemes such as Frandsen [14] because the variables are constructed to automatically satisfy Laplace's equation, and hence the computation is reduced to the solution of two explicit PDEs for the free-surface evolution, making the scheme computationally fast. An essential part of this scheme is the use of the Hilbert-Garrick transform which provides a mapping between the real and imaginary parts of the boundary values of both the mapping and the complex potential for the fluid. Numerical results obtained with this approach are compared to those of Frandsen [14] and they show excellent agreement-when the conformal modulus is properly included in the calculation. This time-dependent conformal mapping strategy is extended to simulate sloshing in a vessel undergoing pendular rotation in Turner et al. [28].

The time-dependent conformal modulus problem
In this section the time-dependence of the conformal modulus is studied for the case where the free surface is specified. First general aspects of time-dependent conformal mapping for domains of the form shown in Fig. 4 are presented. Then two fundamental issues are addressed: the short time behaviour of the conformal modulus, and the implications for error if the conformal modulus is fixed rather than computed. In order to construct the time-dependent conformal mapping which maps the freesurface problem to a rectangle with periodic boundary condition, first form the even extension of the physical problem from x ∈ [0, 1] to x ∈ [0, 2]. The conformal mapping maps the domain x ∈ [0, 2], y ∈ [0, η(x, t)] to μ ∈ [0, 2] and ν = [−Q(t), 0] where Q(t) is the time dependent conformal modulus. Although it is not essential, it will be assumed that the free surface is a graph: y = η(x, t).
In the mapped domain the physical coordinates x and y are harmonic (1.4) and satisfy the Cauchy-Riemann equations (1.3). When these coordinates are evaluated on the free-surface, corresponding to ν = 0 in the mapped domain, they give a parametric form of the free surface The most general solutions for x(μ, ν, t) and y(μ, ν, t) which satisfy the Cauchy-Riemann equations and satisfy condition 1 above are y(μ, ν, t)=B 0 (t)+ν + ∞ n=1 (B n (t)cosh(nπν)+A n (t)sinh(nπν)) cos(nπμ), (2.5) x(μ, ν, t) = μ + ∞ n=1 (A n (t) cosh(nπν) + B n (t) sinh(nπν)) sin(nπμ), (2.6) where A n (t), B n (t) are arbitrary functions of time. Further simplification is obtained by using condition 2 that y = 0 at ν = −Q, which leads to the set of equations Therefore the functions, satisfy conditions 1 and 2. The undetermined functions A n (t) and Q(t) are found by satisfying condition 3. Taking the limit ν → 0 in Eqs. 2.7 and 2.8, the expressions for the free-surface variables are A n (t) tanh (nπQ(t)) cos(nπμ). (2.9)

Small-time asymptotic behaviour of Q(t)
An elementary example illustrating the time dependence of the conformal modulus is obtained by taking a free-surface profile of the form (2.10) At t = 0 the surface is horizontal, but the fluid moves towards x = 1 for t > 0, with a straight free-surface profile. To keep it realistic, the amount of fluid in the vessel is conserved for all time. Although this profile is a simplification to a true sloshing profile, near straight free-surface profiles have been observed in the pendulum vessel simulations of Turner et al. [28].
Substituting the surface Fourier series (2.9) into Eq. 2.10 shows that the functions A n (t) and Q(t) are required to satisfy A n (t) tanh(nπQ(t)) cos(nπμ) Integrating the first expression with respect to μ from 0 to 1 gives which is a relation between the conformal modulus and the odd Fourier modes (The same expression is derived if you integrate from μ = 0 to 2). A second expression linking the A n 's can be derived by multiplying the free surface condition by cos(nπμ) and integrating with respect to μ from 0 to 2. This leads to the equation The A n s and Q can be computed numerically for an arbitrary value of t by considering N terms in the expansions in Eq. 2.9, and solving Eq. 2.12 with Q replaced by Eq. 2.11 via Newton iterations. Finally Q is recovered via Eq. 2.11. The conformal modulus for this numerical approach with N = 500 and t = 1 was checked against the method on page 82 of Papamichael and Stylianopoulos [25] with l = 1/2 which utilizes elliptic integrals, and the results were found to agree.
However, as we are interested in the early time evolution of the conformal modulus, we consider t 1 and calculate an asymptotic solution for the functions A n and substitute these into Eq. 2.11 to give the small time evolution for Q(t).
Substituting the form of Eq. 2.11 into Eq. 2.12 and expanding for small t gives The form of this expansion suggests an expansion for each A n (t) of the form where each A nk ∈ R for all n, k ∈ N/{0}, and therefore at O(t) the leading order term of A n is On substituting this result into Eq. 2.11 we determine the conformal modulus The leading order correction to From the asymptotic form of Q(t) in Eq. 2.13 we note that to leading order and soQ = 0 for all sufficiently small values of t.
With this example we can start to see the implication of fixing Q(t) to be a constant. When Q(t) is fixed, the Fourier coefficients become constrained via Eq. 2.13. In this case, because we have constructed x and y so that they satisfy Laplace's equation, the values of the A n s cannot be found to satisfy both Eqs. 2.11 and 2.12 unless Q is correctly computed.

Induced error due to fixing the modulus
In this subsection additional implications of fixing Q(t) = constant are considered with the more general free surface profile, (2.14) This free surface is the profile associated with the linear solution for the fundamental anti-symmetric free sloshing mode with unit frequency in a fluid of average depth δ. The parameter is a measure of the height of the disturbance to the still fluid depth. At t = 0 the fluid is released from rest and the fluid sloshes back and forth along the vessel.
We can again use the forms of the conformal mapping on ν = 0 given by Eq. 2.9, and we need to calculate values of A n and Q such that Eq. 2.14 is satisfied. This is achieved via a similar approach to that from Section 2.1. Firstly Eq. 2.9 is substituted into Eq. 2.14 giving (2.15) but unlike for the straight free-surface in Section 2.1 it is not possible to derive a simple system of equations for Q(t) and A n (t) by integrating Eq. 2.15. Instead we solve for these functions by choosing Eq. 2.15 to be correct at the N + 1 points, Discretizing in this way (2.15) leads to the N + 1 nonlinear algebraic equations which can be solved via Newton iterations. Before presenting results for general values of , we first consider small amplitude disturbances, 1 and look at the asymptotic structure of the solution, as this will give us an insight into the properties of Q(t).
For small amplitude disturbances we consider the free-surface expression (2.15) with 1. The form of Eq. 2.15 suggests expansions of the form Substituting these into Eq. 2.15 and equating powers of we find that , (2.16) and so Expressions for the coefficients A ij (t) can also be computed to leading order, but are not needed for the discussion here (see below in Fig. 7 for plots of A 1 (t) to leading order). The principal observation is that Q(t) oscillates about some mean value with twice the frequency of the sloshing fluid, and because the correction is O( 2 ), fixing Q = constant in a dynamic calculation means that the errors in the solution will grow more rapidly for larger disturbance amplitudes. The above result suggests that as the simulations become more nonlinear, fixing the conformal modulus may increase the proliferation of numerical errors. This observation is confirmed in the calculations of the sloshing problem for the full irrotational Euler equations in the next section. Figure 6 shows that the O( 2 ) correction to Q(t), with δ = 1, is in good agreement with the numerical result calculated with N = 500 up to = 0.1, while Fig. 7 shows that at = 0.2 the first two terms in the expansion for A 1 gives good agreement with the numerical result.

Conjugate harmonic functions in finite depth
One of the principal advantages of conformal mappings is that the real and imaginary parts of boundary values of holomorphic functions like x + iy and φ + iψ can be related by Hilbert transforms or their generalization.
In the infinite depth case, the conjugate boundary values are related by the Hilbert transform. If u(μ, ν) + iv(μ, ν) is an analytic function in the periodic half space, 0 ≤ μ < 2 and ν < 0, then the boundary values at ν = 0, denoted by U(μ) + iV (μ) are related by the Hilbert transform where the overline denotes average over 0 ≤ μ ≤ 2. The Hilbert transform for functions on 0 ≤ μ ≤ 2 is defined by where the notation used is consistent with that in Papamichael and Stylianopoulos [25]. By combining real and imaginary parts, these two formulas condense into The Hilbert transform has been widely used to develop fast transforms of boundary values for simulations of water waves in infinite depth (e.g. [6,7,11,23,26,32]).
However, when the depth is finite the boundary values at the free surface are coupled to the boundary values at the bottom. This problem of connecting the conjugate boundary values in this case was first solved by Garrick [16] and a detailed derivation is given in Section 1.8 of Papamichael and Stylianopoulos [25] (see also Section 17.4 of Henrici [19] and the paper of Gaier and Papamichael [15]). The theory there is for an annulus, but it carries over with a few modifications to the periodic strip. A sketch of the theory is given here.
Let u(μ, ν) + iv(μ, ν) be an analytic function in the periodic domain 0 ≤ μ < 2 and − Q < ν < 0 , and introduce boundary values where U 0 (μ) and U b (μ) are considered given. The strategy for determining the conjugate functions starts with a solution of the Dirichlet problem for u(μ, ν) with periodicity in μ and the above prescribed boundary values at ν = 0 and ν = −Q. This solution gives u(μ, ν) in terms of the Fourier coefficients of U 0 and U b . The conjugate function is then obtained using the Cauchy-Riemann equations, up to an additive constant. However, in a doubly connected region, a conjugate function exists only if a solvability condition is satisfied (e.g. Theorem 15.1d in Section 15.1 of Henrici [19]). This solvability condition states that a conjugate harmonic function exists if the conjugate period where is a contour which runs from some ν = ν 0 at μ = 0 to ν = ν 0 at μ = 2 in the conformal mapping domain. This essentially means that, if the conjugate function exists, then it is uniquely determined up to an additive constant, which in this case amounts to (Note, the conjugate functions x(μ, ν, t) and y(μ, ν, t) in Eqs. 2.7 and 2.8 do not satisfy this condition, but the conjugate functions x(μ, ν, t) − μ and y(μ, ν, t) − ν do, with U 0 = U b = 0 and V 0 = V b = Q(t). This is utilized in Eq. 5.48). With these two conditions, the conjugate functions u(μ, ν) and v(μ, ν) are completely determined and restriction to the boundaries gives the Hilbert-Garrick transform (3.20) where q = e −πQ , and the new operators are defined by Given a Fourier series representation of a boundary value φ, the action of the above operators is Properties of these operators are summarized in Section 1.8 of Papamichael and Stylianopoulos [25]. An important property is that they are all zero when acting on the constant function.
The above operators satisfy the identities which can be verified by direct calculation. These identities can be used to prove that the inverse map associated with Eq. 3.20 is (3.21) The generalization of the formula (3.19) is In general the boundary values at ν = 0 and ν = −Q are intrinsically linked. The result (3.21) can be used to solve the following problem. Suppose the analytic function z(μ, ν) = u(μ, ν) + iv(μ, ν) is given in the interior with its imaginary part specified on the boundary Im(z) ν=0 = V 0 (μ) and Im(z) with V 0 and V b given. Note, we can neglect the μ-independent contribution from V 0 = V b in Eq. 3.22 as they are not required in Eq. 3.21. Then the real part on each of the boundaries is given by the formula (3.21). Similarly, if the real part is specified on the boundary then the imaginary part is given by the formula (3.20).

The flat bottom case
In general, the real and imaginary parts on the two boundaries are essentially linked. However, in special cases the formulas simplify. For the application in this paper, where the bottom is flat, the imaginary part of the bottom boundary function is constant or zero. Hence the second formula in Eq. 3.20 simplifies to Substitution into the first formula of Eq. 3.20 then gives (3.23) Hence the action of T q on the Fourier series φ is (3.24) The formula (3.24) and its inverse are the key formulas needed for the numerics in the flat bottom case. The reduced operator T q agrees with the operator used in Dyachenko et al. [12], Choi and Camassa [8], Li et al. [22], and Shamin [26]. However, for any more general bottom boundary variation, the full Hilbert-Garrick mappings (3.20) and (3.21) need to be used.

Dynamic sloshing and conformal mapping
Consider dynamic free sloshing in a rectangular vessel for an inviscid, incompressible, irrotational fluid, governed by the Euler equations in the velocity potential formulation. The governing equations for the fluid motion in the interior of the vessel in Fig. 4a can be given in terms of a velocity potential φ and the parametric representation of the free surface X(μ, t) + iY (μ, t). The velocity potential is harmonic φ xx + φ yy = 0 in 0 ≤ x ≤ 2, y ∈ D(t), (4.25) where the fluid domain D(t) is bounded by x = 0, x = 1, y = 0 and above by the unknown surface represented by X(μ, t) + iY (μ, t). The free surface boundary conditions are the dynamic condition where Be(t) is the Bernoulli function, δ is the average fluid depth and we have made the even extension of φ(x, y, t) to x ∈ [1, 2] as in Section 2. The kinematic condition is ∇φ · n = X t · n, on x = X, y = Y, (4.27) where X = (X, Y ) and n is the unit normal at the free surface directed out of the fluid The parametric representation is assumed to be regular: J = 0. The wall boundary conditions are φ y = 0, on y = 0, (4.29) φ x = 0, on x = 0, 1, 2 . (4.30) The moving domain is transformed to a time-dependent rectangle using a conformal mapping, as in Fig. 4b. The mapping x(μ, ν, t) + iy(μ, ν, t) and the complex potential φ(μ, ν, t) + iψ(μ, ν, t) are assumed to be time-dependent holomorphic functions.
We represent the surface variables for φ and ψ with upper case letters Noting that √ J ∇φ · n = − μ in transformed coordinates (cf. equation (A-7) of Bridges and Donaldson [3]), the kinematic boundary condition can be expressed in parametric form as In the transformed domain, the dynamic boundary condition for 0 ≤ μ ≤ 2 is given by after absorbing the Bernoulli function into the potential (cf. equation (3.8) of [3] using Eq. 4.31 above). The bottom boundary condition in transformed coordinates is which has the equivalent representation The free-surface boundary conditions (4.31) and (4.32) are implicit equations for the dependent variables X, Y, and , which makes them unsuitable for standard numerical methods. However, it is possible to reformulate these two equations into three explicit partial differential equations.

Explicit form of the governing equations
For each fixed t, (X(μ, t), Y (μ, t)) is a regular parameterised curve in the plane. Hence any vector in R 2 can be uniquely expressed in terms of the vector (X μ , Y μ ) and its orthogonal vector (−Y μ , X μ ). Express X t in terms of these vectors There are several interesting features to this system. First it is Hamiltonian with the total energy as the Hamiltonian function. The Hamiltonian structure of equations like (4.36)-(4.38) is proved in Benjamin and Bridges [2] for a more general class of two layer flows with differing density. However, the symplectic form is non-constant and nonlinear and so there is no obvious way to devise an effective symplectic numerical integrator (the theory of symplectic integrators relies predominantly on a constant symplectic structure [21]). The second feature of interest is that the α terms are proportional to the tangent vector. Define Z = (X, Y, ), then the explicit equations (4.36)-(4.38) can be written in the compact form The third issue of interest with the form (4.36)-(4.38) is that α provides an automatic time-dependent reparameterization. Introduce an explicit reparameterization Hence taking γ t = α/J eliminates the α term through reparameterization. This implicit reparameterization can be a help or a hindrance. For example this type of reparameterization has been used to enhance the numerical simulation of timedependent free surfaces in Hou et al. [20] and it has been used to enhance the rigorous analysis of the initial value problem for vortex sheets [1].

Using α to enforce analyticity
There is an additional subtlety that arises when Z(μ, t) = X(μ, t) + iY (μ, t) is required to be the boundary value of an analytic function. Let z(ζ, t), ζ = μ + iν, be the mapping function in the interior. Then the boundary values of z t /z ζ are (Note, the right-hand side of the second condition actually yields dQ dt , but this μindependent term term can be neglected as it amounts to the neglected V b term, as discussed after (3.22).) Hence the real part of z t /z μ on each boundary is determined by the Hilbert-Garrick transformation (3.21). Writing (4.34) in complex notation, Z t = ( α + i β)Z μ , β = − μ /J , α = α/J , and using Eq. 3.21, gives α = α + K + R q β , where we write α to mean the mean value of α. Using the definition of T q in Eq. 3.23 the first equation is α = α + T −1 q β , and so the evolution equation for Z is The second equation in Eq. 4.39 is not needed, but provides an interesting identity between the real part at the bottom and the imaginary part at the top.
The function α is not in general arbitrary, it is related to the time-dependent mean part of X(μ, t), as shown by Choi and Camassa [8]. This time-dependent part of X(μ, t) would consist of a function x 0 (t) added to Eq. 2.8. Choi and Camassa [8] then note that x 0 and α, in our notation, are related to one another by taking the mean of the first equation of Eq. 4.40, and that either of these functions can be considered arbitrary, hence implying the other. Taking the mean of the first equation of Eq. 4.40 leads toẋ for a flat bottomed vessel, see Appendix A for workings. Thus, as we have set x 0 = 0 in Eq. 2.8 this implies α = 0. Therefore, the forms of Eqs. 4.36-4.38 we solve are (4.44)
(4.45) Therefore, both φ and ψ can be written as which on the free-surface are given by  Note that in Eqs. 4.42-4.44 ψ only appears as μ , so without loss of generality we can set the function γ 0 (t) to zero.

Numerical scheme
There are several steps in the construction of an algorithm to integrate the initialvalue problem (4.42)-(4.44). Collocation is used to discretize the interval μ ∈ [0, 2] using the 2N points By discretizing using 2N collocation points N terms are retained in each of the infinite summations in Eqs. 2.9 and 4.46. The first derivative terms for the variables in Eqs. 4.42-4.44 are evaluated using second-order, central finite differences. The nonlinearities are treated using pseudo-spectral discretization. The free surface is updated by integrating Eqs. 4.43 and 4.44 using a fourth order Runge-Kutta scheme with a time step t, with specified initial conditions The initial values of X 0 (μ) and Y 0 (μ) at the collocation points along with the initial value of Q come from solving the expression where H (x) is a given free surface elevation, via the numerical scheme described in Section 4. We take advantage of the fact that X(μ, t) and Y (μ, t) are conjugate functions and evolve only the Y t and t equations in Eqs. 4.43 and 4.44 and then X(μ, t) is determined by Similarly does not need to be computed and is determined by where is set to zero via the note at the end of Section 4.3.

Time integration
Time integration of Eqs. 4.43 and 4.44 from t = t m to t = t m+1 proceeds as follows for each time step. Firstly the function T −1 q ( μ /J ) is calculated at the previous timestep, t = t m . Then Eqs. 4.43 and 4.44 are integrated via the Runge-Kutta method to the next time step t = t m+1 . At each internal Runge-Kutta step Y and are calculated via (4.43) and (4.44) and an initial form for X = μ + T −1 q (Y ) and = T q ( ) are determined using Eqs. 2.9 and 4.46 with an initial value of Q given by Q at the previous internal Runge-Kutta step. One could determine X directly by integrating Eq. 4.42, but the numerical scheme is sensitive to numerical errors which occur via this approach. For small t both approaches are equivalent, but for moderate values of t we find integrating Eqs. 4.43 and 4.44 and then finding X via X = μ+T −1 q (Y ) to be more stable, without loss of accuracy.
After updating all four free-surface variables at each internal Runge-Kutta step, we finally update the conformal modulus Q(t). This needs to be done iteratively as the procedure to determine the forms of X(μ, t) and (μ, t) via T q (·) depends upon Q(t). The conformal modulus Q(t) is updated via where B 0 (t) is the non-periodic coefficient of the Fourier transform of Y (μ, t) (see Eq. 2.5) and (n) is the number of iterations. The reason iterations are required is because the correct function X(μ, t) needs to be determined (i.e. the correct value of Q(t)) such that mass conservation in the vessel is achieved. Once Q (n) (t) is updated, X(μ, t) and (μ, t) are updated and this process continues until the relative error in Q (n) (t) is less than some threshold value, which we take to be 10 −10 in Section 6.
Once Q(t) is correctly determined at this internal Runge-Kutta step then the process repeats at the next step. On completion of the composite step of the Runge-Kutta algorithm, Q(t m+1 ) is found via Eq. 5.50, however this final step is not iterative. In order to obtain Q (n) (t) to within a relative error of 10 −10 at each Runge-Kutta step, takes between 3 to 6 iterations on average. The results in Section 6 are presented with N = 400 and t = 1 × 10 −4 . The nonlinear terms are de-aliased, and we apply a filter to X, Y, and after each complete time step to suppress growing higher order Fourier modes not removed by the de-aliasing. These additional higher order modes occur due to the highly nonlinear form of Eqs. 4.43 and 4.44. The reason we are not able to fully de-alias the nonlinear terms, is because the equations contain multiplications by the inverse of a finite Fourier expansion, which itself has an infinite Fourier expansion, and so dealiasing cannot occur completely. However, testing the numerical scheme for various values of N and different de-aliasing values we found that de-aliasing each quadratic nonlinearity using the 2N/3 approach, and applying the additional filtering, to be sufficient to produce converged results.
All results are found to be independent of larger N and smaller t. This can be seen in Fig. 8 which shows values of for the nonlinear result presented in Fig. 10a, where Y (t) = −1 (Y (0, t) − δ). These results show that as t is reduced and N is increased, the results converge and the resulting free-surface plots for N = 300, 400 and 500 are indistinguishable to graphical accuracy (not shown) for 5 × 10 −5 ≤ t ≤ 10 −4 . The slight drop in the value of I for small t values is due to the filtering applied in the numerical scheme to dealias the nonlinear terms. This drop is tolerable for the results presented in this paper, and as stated above, the results are indistinguishable to graphical accuracy.

Numerical results
Here we present linear and nonlinear results from the numerical scheme and compare the linear result with the asymptotics of Section 2.2 and the nonlinear results with the numerical simulations of Frandsen [14]. Frandsen's numerical scheme uses a transformation to map the physical domain to a rectangle for computational ease. The solution for φ is found in the whole domain via an iterative technique, not solely along the free-surface as in our approach. This makes the approach of Frandsen [14] less appealing due to its higher computational cost, but it may be effective in 3D where conformal mapping is no longer viable. As Laplace's equation is solved in the whole domain the resolution of Frandsen's scheme (MM × NN = streamwise points × normal points) cannot be increased to the same resolution as our code without a significant increase in computational time. Hence when using the results of these simulations we shall consider the resolutions 40 × 80 and 80 × 80, both with a time step of 10 −3 , to show that the results have the same qualitative trend for different resolutions. The numerical approach of Frandsen is clearly presented in Frandsen [14] and the reader is directed there for more information. In Frandsen [14], the results of the scheme are compared with third-order asymptotic expansions (in the initial amplitude value) in order to justify the accuracy of the simulations. Thus we have high confidence in the accuracy of this scheme.
The initial conditions for the simulation are taken to be η(x, 0) = δ + cos(π x), and φ(x, y, 0) = 0, identical results to visual accuracy. However, for larger values of this is not the case.
In Fig. 10a we plot Y (t) for = 0.1 for the case when Q is time dependent and the case when we fix Q = 0.4983 which is its initial value. The evolution of Q(t) is given in Fig. 10b. It is clear from Fig. 10a that fixing Q rather than allowing it to be time dependent leads to a different free-surface elevation. This is because the numerical scheme forces the conjugate functions to satisfy Laplace's equation, hence if Q is not correct then an incorrect free-surface is given. For t < 20 the two results are in reasonable agreement, but for larger times the effect of fixing the conformal modulus becomes clearer and more pronounced with a slow drift in the phase of Y (t) and differences in the height of the peaks. However, how can we determine that the result with a time dependent Q is in fact the correct solution? To answer this question we make a comparison with the independent numerical scheme of Frandsen [14] in the next section.  [14] To compare with Frandsen [14], we define the error between the two simulations as the area between the two curves up to t = T , which is given by

Comparison with Frandsen
where η F (x, t) is the free-surface elevation given by Frandsen's numerical scheme.
The motivation behind the simulations in this section are two-fold. Firstly the comparison with an independent numerical scheme will confirm that a time-dependent conformal modulus is necessary in the conformal mapping scheme presented, and secondly it will show the magnitude of the error encountered by fixing Q(t) at the outset of the simulations. In Fig. 11a, b we plot E(45) for various fixed values of the conformal modulus Q fix . We see that E(45) varies as Q fix varies, and has a minimum value at Q fix − δ ≈ −0.005. However, as can be seen in Fig. 11b this minimum value of E(45) is still larger than the error found when Q is time dependent, which is given by the horizontal dashed line. In panels (c) and (d) we plot the same result for the higher resolution run for Frandsen's simulation of 80 × 80 to show that the conclusion is independent of the resolution.
The results in Figs. 10 and 11 suggest that fixing the conformal modulus in dynamic simulations can be an acceptable approximation to the true result, but it depends upon the length of time the simulation is run for, and the level of nonlinearity in the simulation. However, picking the value of Q to minimise this error is not obvious, see Fig. 11, as the obvious choice Q fix = δ does not give the minimum error in this case. In fact, without computing the result with Q = Q(t) it is impossible to know just how good the approximation with Q = Q fix is.

Conclusions
This paper examines the time evolution of the conformal modulus, Q(t), for both a kinematic and dynamic unsteady sloshing problem. By considering the small time evolution of a simplified free-surface between x = 0 and x = 1 modelled as the straight line y = 1 + t (x − 1/2), it was shown that as the profile moves away from a horizontal profile at y = 1, the conformal modulus evolved on an O(t 2 ) timescale. ThusQ = 0 for all times, and therefore it is concluded that performing timedependent dynamic simulations for an inviscid, incompressible, irrotational flow with a fixed conformal modulus will either, lead to results with an incorrect free-surface profile, or lead to a flow which is no longer incompressible and irrotational.
We also presented results of dynamic free sloshing in a rectangular vessel via a numerical scheme which uses time-dependent conformal mappings, and in particular the Hilbert-Garrick transformation. The results of these simulations were compared with simulations using the numerical scheme given in Frandsen [14], and it was found that the minimum error between the two methods occurred when Q was time dependent. 0 = − βY μ + T −1 q ( β)X μ + α, where we have added x 0 (t) to Eq. 2.8. Taking the Fourier representations of X and β, as in Sections 2 and 4, we have where the coefficients are functions of time. Using the formulas given in Section 3 we have the corresponding conjugate functions tanh(kπ Q)(−a k cos(kπ μ)) , coth(kπ Q)(b k sin(kπ μ)) .
Thus using a symbolic algebra package such as MAPLE, it is now easy to show that − βY μ + T −1 q ( β)X μ = 0 , hence Eq. 4.41.