Generalized Helmholtz-Kirchhoff model for two dimensional distributed vortex motion

# Generalized Helmholtz-Kirchhoff model for two dimensional distributed vortex motion

## Abstract

The two-dimensional Navier-Stokes equations are rewritten as a system of coupled nonlinear ordinary differential equations. These equations describe the evolution of the moments of an expansion of the vorticity with respect to Hermite functions and of the centers of vorticity concentrations. We prove the convergence of this expansion and show that in the zero viscosity and zero core size limit we formally recover the Helmholtz-Kirchhoff model for the evolution of point-vortices. The present expansion systematically incorporates the effects of both viscosity and finite vortex core size. We also show that a low-order truncation of our expansion leads to the representation of the flow as a system of interacting Gaussian (i.e. Oseen) vortices which previous experimental work has shown to be an accurate approximation to many important physical flows [9].

123

## 1 Introduction

In this paper we represent solutions of the two-dimensional Navier-Stokes equations as a system of interacting vortices. This expansion, which generalizes the Helmholtz-Kirchhoff model of interacting point vortices in an inviscid fluid, systematically incorporates the effects of both vorticity and finite vortex core size. Furthermore, we give conditions which guarantee the convergence of our expansion. Incompressible viscous flow has two standard analytic representations: a formulation in terms of the primitive velocity and pressure variables, and a formulation in terms of the velocity and vorticity variables [7]. The velocity-vorticity representation has particular advantages when boundaries are unimportant, since vorticity cannot be created or destroyed in the interior of a fluid. The vorticity field can also be directly related to physically observed flow structures such as line and ring vortices.

In two space dimensions, the vorticity field has the additional advantage of reducing to a scalar. An early representation of two-dimensional flow in terms of moving point vortices was developed by Helmholtz-Kirchhoff [4] and by Helmholtz [15] . The point vortex model has been studied extensively - a thorough review of the model and recent developments are described in the monograph by Newton [11]. While the Helmholtz-Kirchhoff point vortex model captures many of the basic physical phenomena observed in two-dimensional rotational flows, experiments with even relatively simple vortex configurations exhibit complications far beyond the point vortex predictions [9]. Additionally, the classical point vortex model neglects the effects of viscosity. However, these experiments also reinforce the idea that in many circumstances the fluid flow may be well approximated by a collection of interacting vortices - albeit vortices with finite core size, subject to the effects of viscosity. A few recent studies of the interaction of viscous vortices can be found in [1, 5, 6, 8, 9, 14]. The main focus of these papers is the merger of two like signed vortices. In [8] the authors use a spatial moment model for 2-D Euler equations which later incorporates weak Newtonian viscosity to derive equations of motions for two like signed vortices. A metastable state is found before merger which consists of two rotating, near-circular, vortices. More recently, it has been conjectured that in a two-vortex system the profiles relax to a pair of gaussian vortices before merging [6] . Thus, it is of interest to extend and generalize the Helmholtz-Kirchhoff point vortex model to a model that incorporates non-zero vortex core size and viscous effects while retaining its basic form. Such an extension is the goal of this paper.

The governing equations for the velocity and pressure variables are

 ∂u∂t+(u⋅∇)u=−∇pρ+νΔu, (1)
 ∇⋅u=0, (2)

where is the fluid density and is the kinematic viscosity. Taking the curl of (1) and (2) gives

 ∂ω∂t+(u⋅∇)ω−(ω⋅∇)u=νΔω, (3)
 ω=∇×u,∇⋅ω=0, (4)

which are the governing equations for the velocity-vorticity variables. For two-dimensional flows, the vorticity vector is perpendicular to the plane of the flow, and the third term on the left-hand side of (3) vanishes. The condition is identically satisfied, and equation (3) then reduces to the single scalar equation:

 ∂ω∂t+(u⋅∇)ω=νΔω, (5)

where is the single, non-zero component of the vorticity. A drawback of the formulation in (5) is that the velocity of the fluid is still present in the equation. However, assuming that the vorticity field is sufficiently localized, the velocity vector can be computed in terms of the vorticity by the Biot-Savart law

 u(x)=12π∫R2(x−y)⊥|x−y|2ω(y)dy , (6)

where for a two-vector , .

In this paper we use equations (5) and (6) to develop a vorticity representation of two-dimensional viscous flow. Our representation is based on a decomposition of the vorticity field into a set of moving distributed vortices. Differential equations are derived for the motion of the vortex centers and for the time evolution of the vortex distributions. The evolution of each individual vortex is represented as an expansion with respect to a sequence of Hermite functions. Such expansions have proven useful in theoretical studies of two-dimensional fluid flows ([2], [3]) and the leading order term in this expansion is precisely the Gaussian vortex (i.e. Oseen vortex [13],[12]) whose utility as an approximation for vortex interaction was shown in [9]. We show that the coefficients in this expansion satisfy a system of ordinary differential equations whose coefficients can be explicitly represented in terms of a fixed, computable kernel function. We also prove the convergence of this expansion. It is shown that our representation reduces in the appropriate limit to the Helmholtz-Kirchhoff model, and allows at the same time arbitrarily complex evolution and interaction of the moving vortices.

In the present paper we concentrate on the mathematical formulation of the generalized Helmholtz-Kirchhoff model. In future work we will explore the predictions of this model both numerically and analytically in a number of different physical settings.

## 2 The “multi-vortex” expansion

In this section we separate the solution of the vorticity equation into components and derive separate evolution equations for each component. From a physical point of view this decomposition will be most useful when each of the components corresponds to a localized region of vorticity (e.g. a vortex) well separated from the other lumps but the mathematical development described below is well defined without regard to these physical considerations. However, with this application in mind, we will often refer to each of the components as a “vortex”.

Consider the initial value problem for the two-dimensional vorticity equation

 ∂ω∂t=νΔω−u⋅∇ω, ω=ω(x,t), x∈R2, t>0 (7) ω(x,0)=ω0(x)

where is the velocity field associated to the vorticity field . We begin by decomposing the initial vorticity distribution by writing

 ω0(x)=N∑j=1ωj0(x) . (8)

Of course this decomposition is not unique - even the number of pieces, , into which we decompose the vorticity is up to us to choose. In general the choice we make will be motivated by physical considerations, however, for the development below, all we require of the decomposition is that the total vorticity of each vortex is non-zero, i.e.

 mj=∫R2ωj0(x)dx≠0, j=1,…,N. (9)

If (9) is satisfied we define by

 ∫R2(x−xj0)ωj0(x)dx=0 , (10)

or equivalently

 xj0=1mj∫R2xωj0(x)dx . (11)

We now write the vorticity for as

 ω(x,t)=N∑j=1ωj(x−xj(t);t) (12)

and the velocity field as

 u(x,t)=N∑j=1uj(x−xj(t);t) (13)

where is the velocity field associated to by the Biot-Savart Law. Of course we still have to define the equations of motion for and .

The centers of the vorticity regions, , and the vorticity regions themselves evolve via a coupled system of ordinary-partial differential equations constructed so that in the limit of zero viscosity and when the different components of the vorticity happen to be point vortices (i.e. Dirac-delta functions) we recover the Helmholtz-Kirchhoff point vortex equations. If we take the partial derivative of (12) and use the equation satisfied by the vorticity, we find:

 ∂tω(x,t) = N∑j=1∂tωj(s−sj(t),t)−N∑j=1˙xj(t)⋅∇ωj(x−xj(t),t) = N∑j=1νΔωj(x−xj(t),t)−N∑j=1(N∑ℓ=1uℓ(x−xℓ(t),t))⋅∇ωj(x−xj(t),t) .

Given this equation it is natural to define as the solution of the equation:

 ∂ωj∂t(x−xj(t),t) = νΔωj(x−xj(t),t)−(N∑ℓ=1uℓ(x−xℓ(t),t))⋅∇ωj(x−xj(t),t) (15) +˙xj(t)⋅∇ωj(x−xj(t),t) , j=1,…,N.

To close this system of equations we must specify how the centers of vorticity evolve. We impose the condition that the first moment of each vorticity region must vanish at every time , i.e. we require that

 ∫R2(x−xj(t))ωj(x−xj(t),t)dx=0  for all  t>0, j=1,…N . (16)

(Note that this equation really contains two conditions - one for each component of .) We impose this condition to fix the evolution of because Gallay and Wayne have recently shown [2] that if one considers the evolution of general solutions of (2) the solution will approach an Oseen vortex, and the rate of the approach will be faster if the vorticity distribution has first moment equal to zero. Solutions of (2) preserve the first moment, and hence if the initial conditions have first moment equal to zero the solution will have first moment zero for all time. The equations (15) no longer preserve the first moment and thus we impose this condition for all time, which then defines the motion of the center of vorticity.

Note that if we change variables in (16) to , we find

 ∫R2zωj(z,t)dz=0 . (17)

Since this equation holds for all we can differentiate both sides with respect to to obtain

 ∫R2z∂tωj(z,t)dz=0 . (18)

Using (15) we can insert the formula for into this integral and we obtain:

 +∫R2z(˙xj(t)⋅∇ωj(z,t))dz=0 . (19)

We first note that if we integrate twice by parts, we have

 ∫R2zΔωj(z,t)dz=0 . (20)

Next if we take the component and integrate by parts we find

 ∫R2z(˙xjn(t)⋅∇ωj(z,t))dz=˙xjn(t)∫R2z∂znωj(z,t)dz=−mj˙xjn(t) , (21)

where and .

###### Remark 2.1.

Note that the equations (15) do preserve the total integral (“mass”) of the solution so this definition of is consistent with (9).

Finally, recalling that the velocity field is incompressible we can rewrite

 (N∑ℓ=1uℓ(z+xj(t)−xℓ(t),t))⋅∇ωj(z,t)=∇⋅(N∑ℓ=1uℓ(z+xj(t)−xℓ(t),t)ωj(z,t))

so, again considering the component and integrating by parts we have

 ∫R2zn(N∑ℓ=1uℓ(z+xj(t)−xℓ(t),t))⋅∇ωj(z,t)dz (22) =−∫R2(N∑ℓ=1uℓn(z+xj(t)−xℓ(t),t))ωj(z,t)dz

Thus, if we combine (20) (21), and (22) we see that (2) reduces to the system of ordinary differential equations for the centers of the vorticity distributions:

 dxjndt(t)=1mjN∑ℓ=1∫R2(uℓn(z+xj(t)−xℓ(t),t)ωj(z,t))dz , (23)

supplemented by the initial conditions (10), while the components of the vorticity evolve according to the partial differential equations (15) with initial conditions

 ωj(z,0)=ωj0(z+xj0) (24)

obtained by combining (8) and (12).

###### Remark 2.2.

Consider (23) in the limit in which the components are all point vortices, i.e. , with the Dirac-delta function. Recall that the velocity field associated with such a point vortex is

 Un(z1,z2,t)=−2∑j=1ϵn,jzj1(z21+z22)

where is the antisymmetric tensor with two indices. Then if we ignore the (singular) term with in the sum on the right hand side of (23) we find

 dxjndt(t)=N∑ℓ=1;ℓ≠jmℓ(∑2k=1ϵn,k(xj(t)−xℓ(t))k)|xj(t)−xℓ(t)|2 (25)

These of course are just Helmholtz-Kirchhoff equations for the inviscid motion of a system of point vortices. Thus, our expansion can be regarded as a generalization of the this approximation which allows for both nonzero viscosity and vortices of finite size. To justify omitting the term with on the right hand side of (25) we note that if we approximate the delta function with a narrow, Gaussian vorticity distribution, and by the corresponding velocity field, this term will vanish by symmetry.

## 3 The moment expansion; case of a single center

In this section we introduce another idea - an expansion of the vorticity in terms of Hermite functions. Then, in the next section we will combine the Hermite expansion with the multi-vortex expansion of the previous section.

The moment expansion is an expansion of the solution of the vorticity equation in terms of Hermite functions. Define

 ϕ00(x,t;λ)=1πλ2e−|x|2/λ2 (26)

where . Three simple facts that we will use repeatedly are

1. Finally, and crucially for what follows, the vorticity function is an exact solution (called the Oseen, or Lamb, vortex) of the two dimensional vorticity equation for all values of .

Note that we will often supress the dependence of on when there is no fear of confusion.

We now define the Hermite functions of order by

 ϕk1,k2(x,t;λ)=Dk1x1Dk2x2ϕ00(x,t;λ) (27)

and the corresponding moment expansion of a function by

 ω(x,t)=∞∑k1,k2=1M[k1,k2;t]ϕk1,k2(x,t;λ) (28)

Note that if the function in (28) is the vorticity field of some fluid the linearity of the Biot-Savart law implies that we can expand the associated velocity field as:

 V(x,t)=∞∑k1,k2=1M[k1,k2;t]Vk1,k2(x,t;λ) (29)

where

 Vk1,k2(x,t;λ)=Dk1x1Dk2x2V00(x,t;λ) (30)

and is the velocity field associated with the Gaussian vorticity distribution - explicitly we have:

 V00(x,t;λ)=12π(−x2,x1)|x|2(1−e−|x|2/λ2) , (31)

We define the Hermite polynomials via their generating function:

 Hn1,n2(z;λ)=((Dn1t1Dn2t2e(2t⋅z−t2λ2))|t=0 (32)

Note that the “standard” Hermite polynomials correspond to taking .

Then using the standard orthogonality relationship for the Hermite polynomials:

 ∫R2Hn1,n2(z;λ=1)Hm1,m2(z;λ=1)e−z2dz=π2n1+n2(n1!)(n2!)δn1,m1δn2,m2 (33)

we see that the coefficients in the expansion (28) are defined by the projection operators:

 M[k1,k2;t]=(Pk1,k2ω)(t)=(−1)k1+k2λ2(k1+k2)2k1+k2(k1!)(k2!)∫R2Hk1,k2(z;λ)ω(z,t)dz (34)

### 3.1 Convergence of the moment expansion

In this subsection we derive a criterion for the convergence of the moment expansion derived above and we show that if this criterion is satisfied for then it is satisfied for all subsequent times .

Our convergence criterion is based on the observation that the Hermite functions are, for any value of , the eigenfunctions of the linear operator

 Lλψ=14λ2Δψ+12∇⋅(xψ) . (35)

This fact can be verified by direct computation and is related to the fact that can be transformed into the Hamiltonian quantum mechanical harmonic oscillator.

The Gaussian function plays a crucial role in the convergence proof, and its dependence on the parameter is particularly important in this discussion, so for this subsection only we will define

 ϕ00(x,t;λ)=Φλ(x,t)

to empasize this dependence.

If one now proceeds as in Lemma 4.7 of [3] one can prove

###### Proposition 3.1.

The operator is self-adjoint in the Hilbert space

 Xλ={f∈L2(R2) | Φ−1/2λf∈L2(R2)}

with innerproduct .

An immediate corollary of this proposition and the general theory of self-adjoint operators is

###### Corollary 3.2.

The eigenfunctions of form a complete orthogonal set in the Hilbert space .

and as a corollary of this result and the observation that the eigenfunctions of are precisely our Hermite functions , we have finally

###### Proposition 3.3.

Suppose that

 ∥f∥2λ=∫R2Φ−1λ(x)|f(x)|2dx<∞ ,

then the expansion

 f(x)=∑k1,k2M[k1,k2]ϕk1,k2(x)

converges with respect to the norm on the Hilbert space .

Thus, the following criterion guarantees that the expansion (28) for the vorticity converges:

 ∫R2Φ−1λ(x)(ω(x,t))2dx<∞ . (36)

The main result of this subsection is the following theorem which proves that if our initial vorticity distribution satisfies (36) for some , then the solution of the vorticity equation with that initial condition will satisfy (36) for all time with and hence as a corollary if the initial vorticity distribution satisfies (36) then our moment expansion converges for all times .

###### Theorem 3.4.

Define the weighted enstrophy function

 E(t)=∫R2Φ−1λ(x)(ω(x,t))2dx .

If the initial vorticity distribution is such that for some , and is bounded (in the norm) then is finite for all times .

Proof: The idea of the proof is to derive a differential inequality for which guarantees that if is finite then will be finite for all . Differentiating we obtain

 dEdt(t)=4νλ2E(t)−4νλ4∫R2|x|2Φ−1λ(x)(ω(x,t))2dx (37) +2∫R2Φ−1λ(x)ω(x,t)∂tω(x,t)dx =4νλ2E(t)−4νλ4∫R2|x|2Φ−1λ(x)(ω(x,t))2dx (38) +2∫R2Φ−1λ(x)ω(x,t)(νΔω−u⋅∇ω)dx

We now consider the last term in (38). First note that upon integration by parts we have:

 2∫R2Φ−1λ(x)ω(x,t)(νΔω(x,t))dx=−2ν∫R2Φ−1λ(x)(|∇ω|2+2λ2ωx⋅∇ω)dx . (39)

The right hand side of the expression in (39) can again be broken up into two pieces and the second can be bounded by:

 2ν∫R2Φ−1λ(x)(2λ2ωx⋅∇ω)dx≤ν∫R2Φ−1λ(x)|∇ω|2dx+4νλ4∫R2Φ−1λ(x2ω2)dx . (40)

Finally, we bound the last term in (37), which comes from the nonlinear term in the vorticity equation. In this estimate we use the fact (see [2], Lemma 2.1) that the norm of the velocity field can be bounded by a constant times the sum of the and norms of the vorticity field - i.e by . This observation, combined with the fact that , which is a consequence of the maximum principle, implies that

 ∥u(⋅,t)∥L∞(R2)≤C(∥ω0∥L1(R2)+∥ω0∥L∞(R2)) .

and hence that

 2∫R2Φ−1λ(x)ω(x,t)(u⋅∇ω)dx≤2C(ω0)∫R2Φ−1λ(x)|ω(x,t)||∇ω|dx (41) ≤4C(ω0)ν∫R2Φ−1λ(x)(ω(x,t))2dx+ν∫R2Φ−1λ(x)|∇ω|2dx

If we now combine the inequalities in (39), (40) and (41), with the expression for in (38) we obtain:

 dEdt(t)≤(4C(ω0)ν+4νλ2)E(t) , (42)

from which we see immediately that if is bounded, remains bounded for all time.

### 3.2 Differential equations for the moments

Assuming that the function is a solution of (2), we can derive differential equations satisfied by the moments in (28). Surprisingly the expressions for the coefficients in these expansions are quite simple and explicit.

If we differentiate (28) and assume that is a solution of the two-dimensional vorticity equation we obtain:

 ∂tω=∞∑k1,k2=1dM[k1,k2;t]dtϕk1,k2(x,t;λ)+∞∑k1,k2=1M[k1,k2;t]∂tϕk1,k2(x,t;λ) =∞∑k1,k2=1M[k1,k2;t](νΔϕk1,k2(x,t;λ)) (43) −⎛⎝∞∑ℓ1,ℓ2=1M[ℓ1,ℓ2;t]Vℓ1,ℓ2(x,t;λ)⎞⎠⋅∇⎛⎝∞∑k1,k2=1M[k1,k2;t]ϕk1,k2(x,t;λ)⎞⎠

From the first of the “simple facts” we stated about , we see that the last term on the first line cancels the middle line and hence if we apply the projection operators defined in (34), we are left with the system of ordinary differential equations for the moments

 dM[k1,k2;t]dt = −Pk1,k2⎡⎣⎛⎝∞∑ℓ1,ℓ2=1M[ℓ1,ℓ2;t]Vℓ1,ℓ2(x,t;λ)⎞⎠ ⋅∇(∞∑m1,m2=1M[m1,m2;t]ϕm1,m2(x,t;λ))]

The somewhat surprising fact which, in our opinion, makes the preceding straightforward calculations interesting is that the projection on the right hand side of (3.2) can be computed explicitly in terms of the derivatives of an relatively simple function.

We now explain how this is done. First recall that:

 ϕm1,m2(x,t;λ)=Dm1x1Dm2x2ϕ00(x,λ), Vℓ1,ℓ2(x,t;λ)=Dℓ1x1Dℓ2x2V00(x,λ) (45)

In order to avoid confusing the two sets of derivatives we will rewrite these formulas as

 ϕm1,m2(x,t;λ) = (Dm1a1Dm2a2ϕ00(x+a,λ))|a=0, (46) Vℓ1,ℓ2(x,t;λ) = (Dℓ1b1Dℓ2b2V00(x+b,λ))|b=0

Inserting these formulas into the right hand side of (3.2) and using the formula for the projection operator in terms of the integration against a Hermite polynomial we obtain

 dMdt[k1,k2,t]=−(−1)(k1+k2)λ2(k1+k2)2k1+k2(k1!)(k2!)∑ℓ1,ℓ2∑m1,m2M[ℓ1,ℓ2,t]M[m1,m2,t] (47) ×∫R2Hk1,k2(x)(Dm1x1Dm2x2V00(x;λ))⋅∇x(Dℓ1x1Dℓ2x2ϕ00(x;λ))dx =−(−1)(k1+k2)λ2(k1+k2)2k1+k2(k1!)(k2!)∑ℓ1,ℓ2∑m1,m2M[ℓ1,ℓ2,t]M[m1,m2,t] ×Dk1t1Dk2t2Dm1b1Dm2b2Dℓ1a1Dℓ2a2∇a⋅⎛⎝∫R2e(−t21−t22+2t1x1+2t2x2)λ2V00(x+b;λ)ϕ00(x+a;λ)dx⎞⎠|t=0,a=0,b=0 .

The last equality in this expression results from rewriting the Hermite polynomial in terms of its generating function.

The last step in deriving the equations for the moments is to evaluate the integral in the last line of (47). The key step in this evaluation is to recall that for these incompressible flows the velocity field can be written in terms of the derivatives of the stream function and that the Laplacian of the stream function is minus the vorticity. Thus, we can write:

 V00(x+b;λ)=−∇∗b(Δb)−1ϕ00(x+b) , (48)

where .

Inserting this into the integral in (47) we find

 ∫R2e(−t21−t22+2t1x1+2t2x2)λ2V00(x+b;λ)ϕ00(x+a;λ)dx (49) =−∇∗b(Δb)−1∫R2e(−t21−t22+2t1x1+2t2x2)λ2ϕ00(x+b;λ)ϕ00(x+a;λ)dx

Now note that all three factors in the integrand are Gaussians and thus the integral can be evaluated explicitly, and we find

 ∫R2e(−t21−t22+2t1x1+2t2x2)λ2ϕ00(x+b;λ)ϕ00(x+a;λ)dx (50) =12πλ2e−12λ2(a21+a22−2a1b1+b21−2a2b2+b22+2a1t1+2b1t1+t21+2a2t2+2b2t2+t22)

We next compute the expression

 −∇∗b(Δb)−112πλ2e−12λ2(a21+a22−2a1b1+b21−2a2b2+b22+2a1t1+2b1t1+t21+2a2t2+2b2t2+t22) (51)

Recall that given a vorticity field , is the associated stream function and the velocity field associated with . Since the inverse Laplacian and derivatives in (51) act only on the -dependent parts of the expression we need to evaluate

 −∇∗b(Δb)−112πλ2e−12λ2(−2a1b1+b21−2a2b2+b22+2b1t1+2b2t2) (52) =−e12λ2((t1−a1)2+(t2−a2)2)∇∗b(Δb)−112πλ2e−12λ2((b1+(t1−a1))2+(b2+(t2−a2))2)

But

 −∇∗b(Δb)−11πλ2e−12λ2((b1+(t1−a1))2+(b2+(t2−a2))2)

is just the velocity field associated with a Gaussian vorticity distribution (i.e. an Oseen vortex) centered at the point which we know explicitly. Hence, the expression on the right hand side of (49) has the explicit representation:

 12πe−12λ2(a21+a22+t21+t22)e12λ2((t1−a1)2+(t2−a2)2) (53) ×(−(b2+(t2−a2),(b1+(t1−a1))((b1+(t1−a1))2+(b2+(t2−a2))2)(1−e−12λ2((b1+(t1−a1))2+(b2+(t2−a2))2))

If we now return to (47) we see that in order to compute the coefficients in the moment equations we need to evaluate the divergence of this last expression with respect to which gives

 K(a1,a2,b1,b2,t1,t2;λ)= −⎛⎜ ⎜⎝(a2t1−b2t1+(−a1+b1)t2)π(a12+a22+b12+b22+2b1t1+t12−2a1(b1+t1)