TThis draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics Planar potential flow on Cartesian grids
Diederik Beckers and Jeff D. Eldredge † Mechanical and Aerospace Engineering, University of California, Los Angeles, CA 90095-1597 USA
Potential flow has many applications, including the modelling of unsteady flows in aerodynamics.For these models to work efficiently, it is best to avoid Biot-Savart interactions between thepotential flow elements. This work presents a grid-based solver for potential flows in twodimensions and its use in a vortex model for simulations of separated aerodynamic flows. Thesolver follows the vortex-in-cell approach and discretizes the streamfunction-vorticity Poissonequation on a staggered Cartesian grid. The lattice Green’s function is used to efficiently solvethe discrete Poisson equation with unbounded boundary conditions. In this work, we use severalkey tools that ensure the method works on arbitrary geometries, with and without sharp edges.The immersed boundary projection method is used to account for bodies in the flow and theresulting body forcing Lagrange multiplier is identified as a discrete version of the bound vortexsheet strength. Sharp edges are treated by decomposing the body-forcing Lagrange multiplierinto a singular and smooth part. To enforce the Kutta condition, the smooth part can then beconstrained to remove the singularity introduced by the sharp edge. The resulting constraints andKelvin’s circulation theorem each add Lagrange multipliers to the overall saddle point system.The accuracy of the solver is demonstrated in several problems, including a flat plate sheddingsingular vortex elements. The method shows excellent agreement with a Biot-Savart method whencomparing the vortex element positions and the force.
1. Introduction
Potential flow plays an important role in aerodynamic modelling, but also appears in otherareas such as the modelling of water waves or wind farms and added mass calculations. Besidesits prominent use for steady flow around airfoils at high Reynolds numbers, potential flow theoryhas long provided the tools for vortex methods to simulate unsteady flows around airfoils andbluff bodies. These vortex methods discretize the vorticity in the flow with singular elementssuch as point vortices, vortex sheets, or a combination of both, and, in the case of an inviscidand incompressible model, the irrotational flow outside of these singular vortex elements is apotential flow. The singular vortex elements representing the free vorticity in an inviscid vortexmodel are tracked as Lagrangian points that are advected by the local flow velocity according toHelmholtz’s second theorem. In case the vortex method inserts new vortex elements in the flowbehind a bluff body or at sharp edges, Kelvin’s circulation theorem dictates that the circulationshould be conserved. Singular potential flow elements also serve to enforce the no-penetrationcondition, with the most common choice in vortex methods being a distribution of singularvorticity on the body, denoted as the bound vortex sheet. Besides their choice for the type ofsingular elements, potential flow solvers for vortex methods differ in their way of calculatingthe flow velocity. This can be done by either using direct interaction between the potential flowelements or by calculating the velocity on a grid over the entire domain and the eventual choicedictates the treatment of boundary and edge conditions.In the first approach, the Green’s function of the Laplacian is applied to the Poisson equationin the velocity-vorticity formulation to give the Biot-Savart integral, which provides the exact † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] F e b solution for a velocity field that satisfies unbounded boundary conditions. Biot-Savart vortexmethods often smooth the Biot-Savart kernel, equivalent to replacing point vortices by vortexblobs (Chorin & Bernard 1973), to suppress Kelvin-Helmholtz instabilities below a certainwavelength resulting from the interactions between closely spaced vortex elements. The solutiongenerally requires O ( 𝑁 ) operations, with 𝑁 the number of vortex elements, to sum the influencesof each discretized vortex element on every other element. With fast multipole methods, it scalesoptimally as O ( 𝑁 ) but with a large prefactor and overhead cost. Inviscid vortex methods of thiskind can straightforwardly use the potential flow tools to enforce the no-penetration, such asconformal mapping or solving the integral equation for a surface singularity distribution throughanalytical inversion, panel discretization, or with Fourier expansions. For a detailed review of thissubject, the reader is referred to Cottet & Koumoutsakos (2000) and Eldredge (2019).A second approach to calculate the flow velocity follows from the discretization of the Poissonequation in the velocity-vorticity formulation or streamfunction-vorticity formulation on a gridover the domain of interest and is called a vortex-in-cell (VIC) approach, first developed byChristiansen (1973). The procedure requires first to regularize the vorticity onto the grid,then to solve the discrete Poisson equation, and finally to interpolate the velocity (or curl ofthe streamfunction) to the Lagrangian vortex elements. It introduces discretization errors butonly requires O ( 𝑀 log 𝑀 ) operations, with 𝑀 the number of grid points, to solve the Poissonequation with current numerical techniques and O ( 𝑁 ) operations to perform the regularizationand interpolation. Similar to the regularized Biot-Savart kernel, the grid spacing together withthe vorticity regularization scheme determine the cut-off wavelength below which the Kelvin-Helmholtz instabilities get suppressed.The early work on VIC methods focused on inviscid vortex dynamics using Fourier-basedPoisson solvers with Dirichlet or periodic boundary conditions (Meng & Thomson 1978; Baker1979; Couët et al. et al. et al. et al. et al. et al. (2015) to use an iterative Brinkman penalizationmethod, which Spietz et al. (2017) extended to three dimensions.The immersed boundary method and Brinkman penalization method both smear out theinfluence of the interface onto nearby grid points. LeVeque & Li (1994) developed the immersedinterface method to overcome this issue and to obtain a higher spatial order of accuracy thanthe immersed boundary method. The premise of this method is to discretize the jump conditionscaused by the interface with finite differences instead of discretizing the Dirac delta function,and the result is a sharp representation of the interface with a second or higher-order accuracy.Marichal et al. (2014) applied the explicit-jump immersed interface method (Wiegmann & Bube2000) in his potential flow method and condensed the influence of the interface in an extra sourceterm in his streamfunction Poisson equation and recognizes that the term is equivalent to a boundvortex sheet strength regularized to the grid. He presents results from the flow over a cylinderand an airfoil, for which they discretized the Kutta condition on the grid with finite-differences.Gillis et al. (2018) extended the method to three dimensions but applies the immersed interfacemethod on the scalar potential instead. He models the potential flow around a sphere and solvesthe Poisson equation using the lattice Green’s function, which automatically takes care of thefar-field boundary conditions.The potential flow method that we propose in this present work is a VIC method on aCartesian grid that uses point vortices to represent the bulk vorticity. We use the lattice Green’sfunction to solve the streamfunction-vorticity Poisson equation with the immersed boundaryprojection method, which was initially developed for viscous flow by Colonius & Taira (2008).The corresponding Lagrange multiplier in the streamfunction Poisson equation is identified as thediscrete analogue of the bound vortex sheet strength. This turns the attention to how the conceptof the discrete vortex sheet strength can be leveraged to enforce the Kutta condition. Instead ofdiscretizing a continuous version of the Kutta condition, we draw inspiration from the analyticaltreatment of the Kutta condition in Biot-Savart methods and eliminate the singularities fromthe discrete bound vortex sheet to ensure the flow velocity is bounded at a sharp edge. This isequivalent to constraining the algebraic system arising from the immersed boundary projectionmethod to make the overall system well-behaved and can therefore be posed as a linear algebraproblem. The result is a saddle point system that we solve with the Schur’s complement method.
2. Methodology
The basic two-dimensional potential flow problem
We consider here a staggered, Cartesian grid with uniform cell size Δ 𝑥 and of infinite extent.The space corresponding to data at cell vertices (nodes) on this grid is denoted by N , and thephysical coordinates of these nodes by x and y . The basic (unbounded) potential flow problem isexpressed as Ls = − w , (2.1)where L is the discrete 5-point Laplacian operator, s ∈ N is the discrete streamfunction, and w ∈ N the discrete vorticity. Following the vortex-in-cell approach (Christiansen 1973), thediscrete vorticity is obtained by regularizing the vorticity from the 𝑁 𝑣 vortex elements onto thecell vertices, w 𝑖, 𝑗 = 𝑁 𝑣 ∑︁ 𝑞 = Δ 𝑥 Γ 𝑣,𝑞 𝑑 (cid:18) x 𝑖 − 𝑋 𝑞 Δ 𝑥 (cid:19) 𝑑 (cid:18) y 𝑗 − 𝑌 𝑞 Δ 𝑥 (cid:19) , (2.2)where 𝑑 is the 𝑀 (cid:48) interpolation kernel (Monaghan 1985), and ( 𝑋 𝑞 , 𝑌 𝑞 ) and Γ 𝑣,𝑞 are the positionand strength of the 𝑞 th vortex element.The discrete velocity field v , whose components lie on the faces of the cells with thecorresponding normals, is computed from s by the discrete curl operation, v = Cs . (2.3)The operator C applies centered differences between the nodes to obtain the velocity componentsat the intermediate centers of cell faces. We denote the space of data that lie on cell faces by F ,so C : N ↦→ F . After the velocity field is computed, the velocities at the positions of any vortex − − . . − − . . 𝑥 𝑦 − − − − − − Δ 𝑥 E rr o r , 𝜖 ( 𝑎 ) ( 𝑏 ) Figure 1. ( a ) Contours of the discrete streamfunction for randomly positioned point vortices ( • ) of randomstrengths, and ( b ) its error ( ◦ ) over the shaded area for different grid spacings. Overlaid is an error (- - - -)that scales as Δ 𝑥 . element in the domain can be computed through interpolation, ( 𝑈 𝑞 , 𝑉 𝑞 ) = ∑︁ 𝑖, 𝑗 v N 𝑖 𝑗 𝑑 (cid:18) x 𝑖 − 𝑋 𝑞 Δ 𝑥 (cid:19) (cid:18) y 𝑗 − 𝑌 𝑞 Δ 𝑥 (cid:19) , (2.4)where v N is v interpolated to the nodes to obtain an overall interpolation scheme that is consistentwith (2.2).The grid differencing operators L and C , and others not yet defined, only comprise differencesbetween grid values and do not scale their results by the grid spacing. Thus, the discrete fieldsrequire re-scaling by this grid spacing to form approximations to their continuous analogs.By convention, we assume v to serve directly as an approximation for the continuous velocityfield. Thus, s is approximately equal to the continuous streamfunction divided by Δ 𝑥 , and w isapproximately the continuous vorticity multiplied by Δ 𝑥 .The particular solution of equation (2.1) can be written down immediately with the help of thelattice Green’s function for L (Katsura & Inawashiro 1971; Cserti 2000; Liska & Colonius 2014).We denote this simply by the inverse operator, s = − L − w . (2.5)It can be shown that, with a suitable truncation of the grid, both L and its inverse are symmetricoperators. However, L is only positive semi-definite, and an additional homogeneous solution, s ∞ ∈ N —for example, corresponding to a uniform flow—can be added to the particular solutionof this equation − L − w + s ∞ . This homogeneous solution allows us to satisfy boundary conditionsat (discrete) infinity. The size of the domain in our simulations is therefore not relevant, as longas it includes the features that are of interest.Figure 1 shows a mesh refinement analysis for a flow consisting of point vortices of randomstrength that are randomly positioned in the lower-left quadrant of the domain. The analysisverifies that the discretization technique of the Poisson equation is second-order accurate in Δ 𝑥 . We compute the error as 𝜖 = (cid:107) 𝜓 ( x , y )/ Δ 𝑥 − s (cid:107) /(cid:107) s (cid:107) , where 𝜓 is the exact solution forthe streamfunction. We only consider the values in the upper right quadrant of the domain toexclude the positions of the point vortices, because the exact singularities at these positions arenot comparable to the regularized, discrete version.2.2. Potential flow with an impenetrable surface
Now, let us suppose we have a rigid impenetrable surface, on which we seek to enforce theno-penetration condition.2.2.1.
Discrete surface and its immersion in the grid
We carry this out at a finite number 𝑁 of discrete surface forcing points; the space of scalardata on these points is denoted by S 𝑁 . In particular, let us define r 𝑥 , r 𝑦 ∈ S 𝑁 as the vectors of 𝑥 and 𝑦 coordinates of the surface points. Each surface point 𝑝 is associated with a small straightsegment of length 𝛿𝑆 𝑝 . Some of the calculations will require information about the local surfaceorientation. For this purpose, we define vectors n 𝑥 , n 𝑦 ∈ S 𝑁 of components of the discrete surfacenormals multiplied by the length of the surface segment 𝛿𝑆 𝑝 . For convenience, let us also defineunit vectors e 𝑝 on this space, equal to 1 at surface point 𝑝 (1 (cid:54) 𝑝 (cid:54) 𝑁 ) and zero at every otherpoint. For example, the 𝑥 coordinate of point 𝑝 is picked out of the vector r 𝑥 by projection ontothe 𝑝 th unit vector: e 𝑇𝑝 r 𝑥 . (2.6)Another vector we will make substantial use of in this paper is ∈ S 𝑁 , a vector of ones on allsurface points.From any vector v ∈ S 𝑁 , we can also form a diagonal 𝑁 × 𝑁 operator D v with the entries ofthe vector along the diagonal. When this operator acts upon another vector w ∈ S 𝑁 , it representsthe Hadamard (i.e., element-by-element) product of the two vectors, D v w = v ◦ w ∈ S 𝑁 . Notethat D v w = D w v , and that D v = v .Surface data are immersed into the grid with the regularization operator R : S 𝑁 ↦→ N . R represents the same operation as applying (2.2) repeatedly to transfer Γ 𝑣,𝑘 / Δ 𝑥 from all vortexelements to the nodes, but operates on arbitrary surface data instead. Grid data are interpolatedonto the surface points with the interpolation operator E : N ↦→ S 𝑁 . E represents the sameoperation as applying (2.4) repeatedly for all grid points but transfers grid data to surface pointsinstead of vortex elements. It should be noted that E can be constructed (and we will assume ithas) so that it is the transpose of the interpolation operator, E = R 𝑇 .2.2.2. The immersed surface potential flow problem
The surface’s motion is specified by a velocity distribution 𝒗 𝑏 , represented discretely bycomponents v 𝑏,𝑥 , v 𝑏,𝑦 ∈ S 𝑁 . The no-penetration condition asserts that the normal componentsof the fluid velocity and this surface velocity must be equal. For rigid bodies, this surface motioncan be alternatively described by a streamfunction, and the no-penetration condition can beimposed equivalently (in two dimensions) by setting the fluid streamfunction equal to that of thesurface s 𝑏 ∈ S 𝑁 , up to a uniform value, s ∈ S 𝑁 . Specifically, translation at velocity ( 𝑈, 𝑉 ) androtation at angular velocity Ω would be described equivalently by velocity components v 𝑏,𝑥 = 𝑈 − Ω r 𝑦 , v 𝑏,𝑦 = 𝑉 + Ω r 𝑥 (2.7)or by a surface streamfunction s 𝑏 = 𝑈 r 𝑦 − 𝑉 r 𝑥 − Ω (cid:16) D r 𝑥 r 𝑥 + D r 𝑦 r 𝑦 (cid:17) . (2.8)Thus, the discrete no-penetration constraint on streamfunction is Es = s 𝑏 − Es ∞ − s . (2.9)For later shorthand, we will denote the difference between the body motion streamfunction andinterpolated uniform flow streamfunction by s (cid:48) 𝑏 ≡ s 𝑏 − Es ∞ . This modified streamfunction simplyconsists of subtracting the components ( 𝑈 ∞ , 𝑉 ∞ ) of the uniform flow from ( 𝑈, 𝑉 ) in (2.8). Theuniform value s is left unspecified and will later serve the role of enforcing a constraint oncirculation. For now, we will suppose that it can be set arbitrarily.The no-penetration constraint is enforced in the basic potential flow problem (2.1) with thehelp of a vector of Lagrange multipliers, f ∈ S 𝑁 , on the surface points. The modified potentialflow problem is thus Ls + R f = − w . (2.10)In fact, by simple comparison with the vorticity on the right-hand side, it is clear that thevector f represents the strength of the discrete bound vortex sheet on the surface. Suppose weconsider the bound vortex sheet 𝛾 ( 𝑠 ) that emerges from the analogous continuous problem on theundiscretized surface, where 𝑠 is the arc-length parameter along the surface. At each point 𝑝 , thediscrete solution f is approximately equal to this continuous solution, multiplied by the length 𝛿𝑆 𝑝 of the small segment surrounding the point: e 𝑇𝑝 f ≈ 𝛾 ( 𝑠 𝑝 ) 𝛿𝑆 𝑝 . (2.11)Thus, the potential flow problem in the presence of the impenetrable surface is (cid:20) L RE (cid:21) (cid:18) s f (cid:19) = (cid:18) − w s (cid:48) 𝑏 − s (cid:19) . (2.12)This problem (2.12) has the structure of a generic saddle-point problem (Benzi et al. Ls ∗ = − w (2.13) S f = s (cid:48) 𝑏 − s − Es ∗ (2.14) s = s ∗ − L − R f , (2.15)where the Schur complement S is S = − EL − R . (2.16)Based on the properties of the matrices comprising S , this operator is symmetric and negativedefinite, and therefore invertible. Its inverse S − , also symmetric, maps a surface distribution ofstreamfunction to a corresponding bound vortex sheet strength.We can describe this algorithm in words: First, solve for the intermediate streamfunction field,associated with vorticity in the fluid, but without regard for the presence of the surface. Second,find the bound vortex sheet whose associated streamfunction cancels the difference betweenthe specified streamfunction on the surface and the intermediate streamfunction evaluated onthe surface. Finally, correct the intermediate streamfunction field for the influence of the boundvortex sheet. A version of the Julia code that implements this algorithm, as well as the algorithmsin the following sections, is available at the Github repository Beckers & Eldredge (2021).We give two examples of the streamfunction with a body present and show the associatedvortex sheet strength. Figure 2 shows a vortex near a circular cylinder and figure 3 shows acircular cylinder that translates horizontally. In both cases, the vortex sheet strength is in verygood agreement with the exact solution from potential flow theory.For later use, we note that the solution of (2.12) can also be written in inverse form usingequation (A 7): (cid:18) s f (cid:19) = (cid:20) L − + L − RS − EL − − L − RS − − S − EL − S − (cid:21) (cid:18) − w s (cid:48) 𝑏 − s (cid:19) . (2.17)The matrix operator in (2.17) is the inverse of the basic saddle-point system. − − − − 𝑥 / 𝑅 𝑐 𝑦 / 𝑅 𝑐 . . . . − . − . − . − . 𝑘 / 𝑁 f 𝑅 𝑐 / ( 𝛿 𝑆 Γ 𝑣 ) ( 𝑎 ) ( 𝑏 ) Figure 2. ( a ) Contours of the discrete streamfunction for a point vortex ( • ) with strength Γ 𝑣 at ( 𝑅 𝑣 , ) near a circular cylinder consisting of 𝑁 points with radius 𝑅 𝑐 and 𝑅 𝑣 / 𝑅𝑐 = /
2, and ( b ) its scaled discretevortex sheet strength ( ◦ ) in function of the normalized surface point index 𝑘 / 𝑁 . 𝑘 = Δ 𝑥 / 𝑅 𝑐 = .
03 and 𝛿𝑆 / Δ 𝑥 = − − − − 𝑥 / 𝑅 𝑦 / 𝑅 . . . . − − 𝑘 / 𝑁 f / ( 𝛿 𝑆 𝑈 ) ( 𝑎 ) ( 𝑏 ) Figure 3. ( a ) Contours of the discrete streamfunction for a horizontally translating circular cylinder withradius 𝑅 , and ( b ) its scaled discrete vortex sheet strength ( ◦ ). Overlaid is the exact continuous solution(- - - -). The simulation is performed with Δ 𝑥 / 𝑅 = .
03 and 𝛿𝑆 / Δ 𝑥 = Non-uniqueness and discrete circulation
In two-dimensional potential flows, there is no unique solution to problem (2.12), since onecan choose any value for the uniform value s and still enforce the no-penetration condition.Equivalently, we can specify any circulation about the body and still enforce this condition. Letus determine the relationship between s and circulation. For later use, let us write this uniformsurface streamfunction as s = 𝑠 , where 𝑠 is a single scalar value. The discrete circulation Γ 𝑏 about the body is given by the sum of the bound vortex sheet data and can be written compactlyas Γ 𝑏 = 𝑇 f . (2.18) . . . . − − . − − . 𝑘 / 𝑁 f / 𝛿 𝑆 ( 𝑎 ) ( 𝑏 ) Figure 4. ( a ) Geometry and ( b ) the scaled bound vortex sheet strength associated with a uniform,unit-strength streamfunction for elliptical cylinders with different aspect ratios ( 𝐴𝑅 ): · · · · · · , 𝐴𝑅 = 𝐴𝑅 =
2; – – –, 𝐴𝑅 =
3; ——, 𝐴𝑅 = ∞ (flat plate). The simulations are performed with Δ 𝑥 / 𝑅 = . 𝛿𝑆 / Δ 𝑥 = The circulation of the vortex sheet in the solution (2.14) is Γ 𝑏 = 𝑇 S − (cid:16) s (cid:48) 𝑏 + EL − w (cid:17) − 𝑠 𝑇 S − . (2.19)The scalar factor 𝑇 S − in this expression is a property of the set of points and their immersioninto the Cartesian grid. Part of this factor, S − , represents the bound vortex sheet strengthassociated with a uniform, unit-strength streamfunction on the surface. This sheet has a particularlyimportant role in some of the discussion to follow, so we will denote its strength by f : f ≡ S − . (2.20)The transpose of f , equal to 𝑇 S − , calculates the circulation of the associated bound vortexsheet when it acts upon a surface streamfunction. Thus, the factor 𝑇 S − is the circulationassociated with a uniform, unit-strength surface streamfunction. We will refer to this as Γ : Γ ≡ 𝑇 S − ≡ f 𝑇 ≡ 𝑇 f . (2.21)Figure 4 demonstrates the distribution of f for elliptical cylinders with different aspect ratios. Fora circular cylinder, f assumes a uniform distribution, which agrees with the fact that tangentialvelocity around a circular cylinder with non-zero circulation is constant at a given radius. Whenthe aspect ratio increases, the distribution gradually shows stronger variations near the edges ofthe major axis and eventually turns singular for the flat plate. This will be discussed in more detailin the next section.The last term in (2.19) illustrates the direct relationship between the scalar value 𝑠 and thebound circulation Γ 𝑏 , and one can interpret 𝑠 as a means of setting the circulation. Thus far,the value 𝑠 has not been fixed. We will use it in the next section to enforce the Kutta condition.However, we will use it here for an immediate purpose. The prescribed surface streamfunction s (cid:48) 𝑏 (given by (2.8), with the uniform flow accounted for) may have some associated boundcirculation, and it is desirable to adjust it by adding or subtracting a uniform value so that it hasnone. Equation (2.19) suggests that this circulation can be removed by setting 𝑠 to 𝑇 S − s (cid:48) 𝑏 / Γ and then subtracting this value (multiplied by the uniform vector ) from s (cid:48) 𝑏 . Overall, this processcan be encapsulated in a circulation removal operator that acts upon the surface streamfunction P Γ ≡ I − 𝑇 Γ S − . (2.22)It is easy to verify that 𝑇 S − P Γ =
0, so that the circulation of any surface streamfunction actedupon by P Γ is indeed zero. It is important to observe, also, that s (cid:48) 𝑏 can be replaced by P Γ s (cid:48) 𝑏 withoutaffecting the nature of the no-penetration condition. We also note that the composite operator S − P Γ is symmetric, just as S − is: S − P Γ = S − − S − 𝑇 S − Γ = S − − f f 𝑇 Γ . (2.23)2.4. The Kutta condition
For surfaces that contain convex edges, the vortex sheet strength assumes a singular behavior inthe vicinity of these edges, with a strength that depends on the interior angle of the edge: sharperedges have more singular behavior. In the discrete representation of the surface, edges are onlyapproximately represented by the sudden disruptions of positions in clusters of adjacent points.The behavior in this discrete form is not quite singular, but the solution of (2.14) nonethelessexhibits a large and rapid change of amplitude.If we seek to eliminate this behavior, we must first have some means of exposing it. In fact,for any discretized surface, the essence of this nearly-singular behavior lies in the vector f , andall other bound vortex sheets associated with the same surface share the same nearly-singularbehavior. Thus, we will use a multiplicative decomposition of the vortex sheet strength: f = f ◦ ˜ f , (2.24)where ◦ is the Hadamard product. This decomposed form isolates the singular behavior into f ,and ˜ f is a relatively smoother vector of surface point data. In the regularization operation on f ,we can absorb f into R , first noting that the Hadamard product can alternatively be written withthe help of a diagonal matrix, f = f ◦ ˜ f = D f ˜ f . (2.25)Then, we can define a re-scaled regularization operator, R f = RD f ˜ f = ˜ R ˜ f . (2.26)The re-scaled operator ˜ R = RD f can, in turn, be absorbed into the Schur complement, defining˜ S = − EL − ˜ R = SD f . A useful property of ˜ S is that it preserves uniform vectors:˜ S = SD f = S f = , ˜ S − = . (2.27)The decomposition of the vortex sheet strength is demonstrated in figure 5 for a flat plate in auniform flow. The original vortex sheet strength shows large amplitude variations at the leadingand trailing edge, corresponding to large variations in the flow velocity field as the flow tries tomake its way around the edges. By use of decomposition (2.24), these discrete singularities areretained in f and we are left with a smooth ˜ f , which varies linearly from the leading edge ( 𝑘 = 𝑘 = 𝑁 ) .The Kutta condition corresponds to annihilating the nearly-singular behavior at a surface point.At such points, we will set the corresponding value of ˜ f to zero. Suppose we wish to enforce theKutta condition at an edge corresponding to surface point 𝑘 . The condition is e 𝑇𝑘 ˜ f = . (2.28)0 − − − f / ( 𝛿 𝑆 𝑈 ∞ ) . . . . − − . . 𝑘 / 𝑁 ˜ f . . . . − − 𝑘 / 𝑁 f / ( 𝛿 𝑆 𝑈 ∞ ) − − . . − − . . 𝑥 / 𝑐 𝑦 / 𝑐 = ×( 𝑎 )( 𝑏 )( 𝑐 ) ( 𝑑 ) Figure 5. ( a ) Contours of the discrete streamfunction for a flat plate of length 𝑐 at 30 ° in a uniform flow 𝑈 ∞ without enforcement of the Kutta condition. The point-wise product of ( b ) the discrete vortex sheet strengthassociated with a uniform, unit-strength streamfunction on the body, and ( c ) the smooth vector that resultsfrom re-scaling the regularization operator, composes ( d ) the discrete vortex sheet strength. The simulationis performed with Δ 𝑥 / 𝑐 = .
01 and 𝛿𝑆 / Δ 𝑥 = Using the Kutta condition in a steady-state problem
We will first take the steady-state approach to enforce the Kutta condition: allow the boundcirculation to be set appropriately, with the implicit understanding that there is a starting vortexof equal and opposite circulation at infinity that preserves the Kelvin circulation theorem. TheLagrange multiplier for this constraint will not be Γ 𝑏 , but 𝑠 . We also use the circulation removaloperator to adjust the imposed surface streamfunction: L ˜ R E e 𝑇𝑘 (cid:169)(cid:173)(cid:171) s ˜ f 𝑠 (cid:170)(cid:174)(cid:172) = (cid:169)(cid:173)(cid:171) − wP Γ s (cid:48) 𝑏 (cid:170)(cid:174)(cid:172) (2.29)This block system, like the earlier one in (2.12), has a saddle point form, and we can reduce itby the same block-LU decomposition to develop a solution algorithm. We will interpret it in thegeneral form (A 3), with the upper left 2 × A , the solution vector x and1constraint force y set, respectively, to x = (cid:18) s ˜ f (cid:19) , y = 𝑠 , (2.30)the remaining operators set to B = (cid:2) e 𝑇𝑘 (cid:3) , B 𝑇 = (cid:20) (cid:21) , C = , (2.31)and the right-hand side vectors set to r = (cid:18) − wP Γ s (cid:48) 𝑏 (cid:19) , r = . (2.32)We note that block A has the original form of the system before the Kutta constraint (2.12),though with the slight modification of a re-scaled regularization operator, and we already have theinverse of A available from (2.17). The solution of this original system forms the intermediatesolution of the full system endowed with the Kutta condition:˜ f ∗ = ˜ S − (cid:16) P Γ s (cid:48) 𝑏 + EL − w (cid:17) , s ∗ = − L − (cid:16) w + ˜ R ˜ f ∗ (cid:17) . (2.33)Then, using the general procedure outlined in appendix A, the solution of the full system (2.29)is easy to develop; its Schur complement is simply S = − . (2.34)Applying the general solution equations, and using the property (2.27) to simplify the resultingoperators, it can be shown that the solution is 𝑠 = e 𝑇𝑘 ˜ f ∗ (2.35) (cid:18) s ˜ f (cid:19) = (cid:18) s ∗ ˜ f ∗ (cid:19) − (cid:20) − L − ˜ R (cid:21) e 𝑇𝑘 ˜ f ∗ . (2.36)The entire solution can be written more compactly as˜ f = P 𝐾𝑘 ˜ S − (cid:16) P Γ s (cid:48) 𝑏 + EL − w (cid:17) , (2.37) s = − L − (cid:16) w + ˜ R ˜ f (cid:17) , (2.38)where we have defined the Kutta projection operator , P 𝐾𝑘 ≡ I − 𝑇𝑘 , (2.39)which acts upon the (smooth part of the) bound vortex sheet vector, subtracting the value at point 𝑘 from every point, including at 𝑘 itself.The application of the Kutta condition to a steady-state problem is demonstrated in Figure 6on the flat plate problem that was introduced in the previous section. By constraining the trailingedge point of ˜ f , its whole distribution is shifted upward such that the last value equals zero. Theresulting streamfunction indicates that the flow then indeed smoothly leaves the trailing edge.Note that the Lagrange multiplier for the Kutta condition takes the simple value given byequation (2.35), revealing that the additional streamfunction on the surface is exactly the value ofthe intermediate bound vortex sheet at the Kutta point 𝑘 .2.4.2. Using the Kutta condition to set a new vortex element
In the previous section, we used the Kutta condition to set the bound circulation but did notexplicitly create a new vortex element. This vortex element was assumed to lie at infinity so thatits effect was negligible except insofar as it left equal but opposite circulation about the body.2 − − . . − − . . 𝑥 / 𝑐 𝑦 / 𝑐 . . . . . . 𝑘 / 𝑁 ˜ f ( 𝑎 ) ( 𝑏 ) Figure 6. ( a ) Contours of the steady, discrete streamfunction for a flat plate of length 𝑐 at 30 ° in a uniformflow with enforcement of the Kutta condition at the trailing edge, and ( b ) the smooth part of the associateddiscrete vortex sheet strength. The simulation is performed with Δ 𝑥 / 𝑐 = .
01 and 𝛿𝑆 / Δ 𝑥 = In this section, we will create a new vortex element in the vicinity of the edge at which we areapplying the Kutta condition. We will thus seek to establish the strength of this new element andto do so in such a manner that the overall circulation of the flow is conserved. Once the elementis created, it will be allowed to advect with the local fluid velocity (minus its own contribution tothis velocity).Let us assume that the new vortex element (which we label with the subscript 1) is introducedat some point in physical space, and that its immersion into the Cartesian grid is described by agrid vector d ∈ N and that its strength (i.e., its circulation) is 𝛿 Γ . Thus, the fluid vorticity afterthis new element’s introduction can be written as w + 𝛿 Γ d . (2.40)The vector d is the discrete Dirac delta function, identical (or similar) to the function thatconstitutes the regularization R and interpolation E matrices.Kutta condition (2.28) is still to be enforced. We also seek to ensure that the total circulationis zero. (We are assuming that the flow has started from rest.) Let us denote the circulation of theexisting fluid vorticity w by Γ w = Δ 𝑥 𝑇 w , (2.41)where ∈ N is a grid vector of ones. Then, the circulation constraint is 𝑇 f + 𝛿 Γ + Γ w = . (2.42)The circulation of the bound vortex sheet f can be re-written in terms of the smooth part of thesheet as 𝑇 f = f 𝑇 ˜ f .With these two constraints, the overall saddle point system of equations is L ˜ R d E e 𝑇𝑘 f 𝑇 (cid:169)(cid:173)(cid:173)(cid:173)(cid:171) s ˜ f 𝑠 𝛿 Γ (cid:170)(cid:174)(cid:174)(cid:174)(cid:172) = (cid:169)(cid:173)(cid:173)(cid:173)(cid:171) − wP Γ s (cid:48) 𝑏 − Γ w (cid:170)(cid:174)(cid:174)(cid:174)(cid:172) . (2.43)Again, the basic saddle-point matrix constitutes the upper left 2 × A and the solution3vector x is as before. The constraint force vector is y = (cid:18) 𝑠 𝛿 Γ (cid:19) , (2.44)and the remaining vectors and operators are now r = (cid:18) − Γ w (cid:19) , B = (cid:20) e 𝑇𝑘 f 𝑇 (cid:21) , B 𝑇 = (cid:20) d (cid:21) , C = − (cid:20) (cid:21) . (2.45)The solution algorithm follows, once again, from the equations in appendix A. After carryingout the block matrix multiplications, it can be shown that the Schur complement (A 2) is the 2 × S = (cid:20) − e 𝑇𝑘 ˜ f − f 𝑇 + f 𝑇 ˜ f , (cid:21) (2.46)where, for convenience, we have defined˜ f = ˜ S − EL − d , (2.47)which represents the (smooth part of the) strength of the vortex sheet that “reacts” to the presenceof a unit-strength vortex d immersed into the grid, canceling that vortex’s induced velocity onthe surface. The term f 𝑇 ˜ f represents this sheet’s bound circulation and e 𝑇𝑘 ˜ f is its contributionto the Kutta condition at point 𝑘 . The problem (A 6) for the constraint forces 𝑠 and 𝛿 Γ is then (cid:20) − e 𝑇𝑘 ˜ f − f 𝑇 + f 𝑇 ˜ f , (cid:21) (cid:18) 𝑠 𝛿 Γ (cid:19) = (cid:18) − e 𝑇𝑘 ˜ f ∗ − Γ w − f 𝑇 ˜ f ∗ (cid:19) . (2.48)The determinant of this Schur complement matrix is − − f 𝑇 P 𝐾𝑘 ˜ f , which represents the negativeof the circulation of the unit vortex and its associated vortex sheet, after the Kutta condition hasbeen enforced on this sheet. It is straightforward then to calculate the strength of the new vortex 𝛿 Γ and the additional uniform surface streamfunction, 𝑠 : 𝛿 Γ = − Γ w − f 𝑇 P 𝐾𝑘 ˜ f ∗ + f 𝑇 P 𝐾𝑘 ˜ f , 𝑠 = e 𝑇𝑘 (cid:16) ˜ f ∗ + 𝛿 Γ ˜ f (cid:17) . (2.49)where the intermediate solution ˜ f ∗ is available from (2.33). From these, we can then obtain thevortex sheet strength and the fluid streamfunction,˜ f = P 𝐾𝑘 (cid:16) ˜ f ∗ + 𝛿 Γ ˜ f (cid:17) , s = − L − (cid:16) w + 𝛿 Γ d + ˜ R ˜ f (cid:17) . (2.50)We now apply the method to the flat plate problem with a point vortex near the trailing edge toenforce the Kutta condition at that edge. We position the point vortex at a distance 3 Δ 𝑡𝑈 ∞ from theedge in the direction of the free stream, perpendicular to the plate. Figure 7 shows that, becauseof the proximity of the point vortex to the flat plate, ˜ f exhibits a quick variation in its value at thesurface points that lie closest to the point vortex. The value at the trailing edge point itself is stillconstrained to zero and the flow again leaves the edge smoothly. This situation corresponds to theflow right after impulsively starting a uniform flow around a flat plate and the point vortex nowrepresents the starting vortex.2.4.3. Applying more than one Kutta condition on a body
Suppose we wish to enforce the Kutta condition at two edges of the body—at points 𝑘 and 𝑘 —instead of one. Each such point has a constraint, e 𝑇𝑘 𝑗 ˜ f = , 𝑗 = , . (2.51)4 − − . . − − . . 𝑥 / 𝑐 𝑦 / 𝑐 . . . . − . . 𝑘 / 𝑁 ˜ f ( 𝑎 ) ( 𝑏 ) Figure 7. ( a ) Contours of the unsteady, discrete streamfunction for a flat plate of length 𝑐 at 30 ° in a uniformflow with release of vorticity into a point vortex ( • ) for enforcement of the Kutta condition at the trailingedge. ( b ) The smooth part of the associated discrete vortex sheet strength. The simulation is performed with Δ 𝑥 / 𝑐 = .
01 and 𝛿𝑆 / Δ 𝑥 = For two such constraints, we need two Lagrange multipliers: the strengths of two new vortices, 𝛿 Γ and 𝛿 Γ , immersed into the grid with d and d , respectively; and we still need the Lagrangemultiplier 𝑠 to ensure that Kelvin’s circulation theorem is also enforced. The system in theprevious section is thus easily generalized to the following: L ˜ R d d E e 𝑇𝑘 e 𝑇𝑘 f 𝑇 (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) s ˜ f 𝑠 𝛿 Γ 𝛿 Γ (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) − wP Γ s (cid:48) 𝑏 − Γ w (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) . (2.52)The system is reduced in the same manner as before, with the same intermediate solutionobtained from the basic system (2.12). Now, the Schur complement problem for the constraintforces takes the form − e 𝑇𝑘 ˜ f e 𝑇𝑘 ˜ f − e 𝑇𝑘 ˜ f e 𝑇𝑘 ˜ f − f 𝑇 + f 𝑇 ˜ f + f 𝑇 ˜ f , (cid:169)(cid:173)(cid:171) 𝑠 𝛿 Γ 𝛿 Γ (cid:170)(cid:174)(cid:172) = (cid:169)(cid:173)(cid:173)(cid:171) − e 𝑇𝑘 ˜ f ∗ − e 𝑇𝑘 ˜ f ∗ − Γ w − f 𝑇 ˜ f ∗ (cid:170)(cid:174)(cid:174)(cid:172) , (2.53)where we have now defined bound vortex sheets associated with each of the two new vortices(with unit strengths): ˜ f 𝑗 = ˜ S − EL − d 𝑗 , (2.54)for 𝑗 = ,
2. It is interesting to note that, if we take the difference between the two Kutta constraints,we obtain (cid:16) e 𝑇𝑘 − e 𝑇𝑘 (cid:17) (cid:16) ˜ f ∗ + 𝛿 Γ ˜ f + 𝛿 Γ ˜ f (cid:17) = . (2.55)It can be shown that this Schur complement problem can be split into (cid:34) + f 𝑇 P 𝐾𝑘 ˜ f + f 𝑇 P 𝐾𝑘 ˜ f + f 𝑇 P 𝐾𝑘 ˜ f + f 𝑇 P 𝐾𝑘 ˜ f (cid:35) (cid:18) 𝛿 Γ 𝛿 Γ (cid:19) = − (cid:32) Γ w + f 𝑇 P 𝐾𝑘 ˜ f ∗ Γ w + f 𝑇 P 𝐾𝑘 ˜ f ∗ (cid:33) (2.56)5 − − . . − − . . 𝑥 / 𝑐 𝑦 / 𝑐 . . . . − . − . . . 𝑘 / 𝑁 ˜ f ( 𝑎 ) ( 𝑏 ) Figure 8. ( a ) Contours of the unsteady, discrete streamfunction for a flat plate of length 𝑐 at 30 ° in a uniformflow with release of vorticity into two point vortices ( • ) for enforcement of the Kutta condition at bothedges. ( b ) The smooth part of the associated discrete vortex sheet strength. The simulation is performedwith Δ 𝑥 / 𝑐 = .
01 and 𝛿𝑆 / Δ 𝑥 = and 𝑠 = (cid:16) e 𝑇𝑘 + e 𝑇𝑘 (cid:17) (cid:16) ˜ f ∗ + 𝛿 Γ ˜ f + 𝛿 Γ ˜ f (cid:17) . (2.57)The latter equation, when combined with (2.55), reveals that the value of the vortex sheet strength˜ f ∗ + 𝛿 Γ ˜ f + 𝛿 Γ ˜ f is the same at both Kutta points and equal to 𝑠 .Equation (2.56) can be solved easily for the strengths of the two new point vortices. Then, thesolution for the vortex sheet strength and streamfunction are˜ f = (cid:16) P 𝐾𝑘 + P 𝐾𝑘 (cid:17) (cid:16) ˜ f ∗ + 𝛿 Γ ˜ f + 𝛿 Γ ˜ f (cid:17) , s = − L − (cid:16) w + 𝛿 Γ d + 𝛿 Γ d + ˜ R ˜ f (cid:17) (2.58)We now apply this method in figure 8 to enforce the Kutta condition at the leading and trailingedge of our flat plate problem. We position a point vortex close to each edge and observe againthat ˜ f shows strong variation at the surface points closest to the two point vortices. The contoursof the streamfunction indicate that the flow indeed leaves the edges smoothly. Like the previouscase, this solution corresponds to the flow right after impulsively starting a uniform flow arounda flat plate, but unlike the previous case, the flow now separates at the leading edge.It should be observed that these solutions are posed in a manner easily extensible to an arbitrarynumber of edges. 2.5. Generalized edge condition
In the previous section, we demonstrated the means of annihilating the (nearly) singularbehavior at edges on a discretized surface. In some cases, our desire is not to annihilate thisbehavior, but simply to keep it within some bounds. In the analytical treatment of potential flowproblems, this objective is served by placing an inequality constraint on the edge suction parameter
Ramesh et al. (2014); Darakananda & Eldredge (2019); Eldredge (2019). That parameter isproportional to the coefficient on the bound vortex sheet strength’s singularity Eldredge (2019),so in this discrete setting, in which we have extracted the singular part of f in the form of f ,we expect the suction parameter to be related to the value of ˜ f at the edge. In fact, by simple6comparison, it can be shown that e 𝑇𝑘 ˜ f = − 𝜋𝑐 Γ 𝜎 𝑘 (2.59)for a flat plate of length 𝑐 , where 𝜎 𝑘 is the suction parameter at the edge corresponding to point 𝑘 . Let 𝜎 min 𝑘 and 𝜎 max 𝑘 denote the minimum and maximum tolerable values of 𝜎 𝑘 at edge 𝑘 . Wethen seek to confine the suction parameter to the range 𝜎 min 𝑘 (cid:54) 𝜎 𝑘 (cid:54) 𝜎 max 𝑘 . This generalized edgeconstraint is placed on the suction parameter of the intermediate sheet ˜ f ∗ . To avoid confusion,we will redefine the bounds based on this smooth part of the vortex sheet rather than 𝜎 𝑘 itself;for this, we define ˜ 𝑓 min 𝑘 = − 𝜋𝜎 max 𝑘 / Γ and ˜ 𝑓 max 𝑘 = − 𝜋𝜎 min 𝑘 / Γ . Thus, we inspect whether thevalue e 𝑇𝑘 ˜ f lies in the range ˜ 𝑓 min 𝑘 (cid:54) e 𝑇𝑘 ˜ f ∗ (cid:54) ˜ 𝑓 max 𝑘 . (2.60)If e 𝑇𝑘 ˜ f lies within this range, then no new vortex is created near the edge (or equivalently, anew vortex of zero strength is created); if e 𝑇𝑘 ˜ f ∗ > ˜ 𝑓 max 𝑘 , then we create a new vortex so that e 𝑇𝑘 ˜ f = ˜ 𝑓 max 𝑘 ; and if e 𝑇𝑘 ˜ f ∗ < ˜ 𝑓 max 𝑘 , then we do the same, but now so that e 𝑇𝑘 ˜ f = ˜ 𝑓 min 𝑘 . Note that theKutta condition simply corresponds to setting ˜ 𝑓 min 𝑘 = ˜ 𝑓 max 𝑘 = e 𝑇𝑘 ˜ f ∗ < ˜ 𝑓 min 𝑘 and e 𝑇𝑘 ˜ f ∗ > ˜ 𝑓 max 𝑘 ;then we solve the system (cid:34) + f 𝑇 P 𝐾𝑘 ˜ f + f 𝑇 P 𝐾𝑘 ˜ f + f 𝑇 P 𝐾𝑘 ˜ f + f 𝑇 P 𝐾𝑘 ˜ f (cid:35) (cid:18) 𝛿 Γ 𝛿 Γ (cid:19) = − (cid:32) Γ w + f 𝑇 P 𝐾𝑘 ˜ f ∗ + Γ ˜ 𝑓 min 𝑘 Γ w + f 𝑇 P 𝐾𝑘 ˜ f ∗ + Γ ˜ 𝑓 max 𝑘 (cid:33) . (2.61)But if, say, ˜ 𝑓 min 𝑘 (cid:54) e 𝑇𝑘 ˜ f ∗ (cid:54) ˜ 𝑓 max 𝑘 , then we set 𝛿 Γ = 𝛿 Γ = − Γ w + f 𝑇 P 𝐾𝑘 ˜ f ∗ + Γ ˜ 𝑓 min 𝑘 + f 𝑇 P 𝐾𝑘 ˜ f . (2.62)The effect of applying these generalized edge conditions to the leading edge of a flat plate isshown in figure 9 for the first instants after impulsively starting a uniform flow. In this and thefollowing simulations, forward Euler is used to advance the system in time. When new pointvortices are inserted, they are placed one-third of the way from the edge to the last released vortexfrom that edge. The positions of the point vortices emanating from the leading edge in figure 9indicate that as 𝜎 max 𝑘 increases, the stream of point vortices is swept back from the edge. At thetrailing edge, the Kutta condition is enforced in each case and the positions of the point vorticesoverlap, as they are not yet influenced by the different situations at the leading edge in these firstinstants. 2.6. Force and the added mass
Impulse-based calculations of force and moment
We will develop a means of calculating the force and moment on the body through the negativerate of change of impulse in the fluid (Eldredge 2019). The continuous expressions for linear and7 − . − . . . . − . − . − . . . . 𝑥 / 𝑐 𝑦 / 𝑐 . . . . − − . . 𝑘 / 𝑁 ˜ f ( 𝑎 ) ( 𝑏 ) Figure 9. Effect of increasing 𝜎 maxLE / 𝑈 ∞ from 0 ( • and ◦ ) to 0 .
15 ( (cid:4) and (cid:3) ) and 0 . (cid:78) and (cid:77) ) on ( a ) thepositions of shedded point vortices and ( b ) the smooth part of the associated discrete vortex sheet strengthfor a flat plate of length 𝑐 at 60 ° , 0 . 𝑈 ∞ . At thetrailing edge, the Kutta condition is enforced. The simulation is performed with Δ 𝑥 / 𝑐 = . 𝛿𝑆 / Δ 𝑥 = Δ 𝑡𝑈 ∞ / 𝑐 = . angular impulse (about the origin) are, in two dimensions, 𝑷 = ∫ 𝑉 𝑓 𝒙 × 𝝎 d 𝑉 + ∫ 𝑆 𝑏 𝒙 × ( 𝒏 × 𝒗 ) d 𝑆 (2.63) 𝚷 = ∫ 𝑉 𝑓 𝒙 × ( 𝒙 × 𝝎 ) d 𝑉 + ∫ 𝑆 𝑏 𝒙 × [ 𝒙 × ( 𝒏 × 𝒗 )] d 𝑆 (2.64)If there is only a single body, then the force and moment (about the origin) exerted by the fluidon that body are given by 𝑭 = − 𝜌 d 𝑷 d 𝑡 , 𝑴 = − 𝜌 d 𝚷 d 𝑡 , (2.65)where 𝜌 is the fluid density. In the two-dimensional applications of this paper, angular impulseand the moment have only a single component, e.g., 𝚷 = Π 𝒆 𝑧 , where 𝒆 𝑧 is the unit vector outof the plane.It should be observed that, by definition, the bound vortex sheet strength 𝛾 is equal to the jumpin tangential velocity between the fluid and the surface, 𝒏 × 𝒗 = 𝛾 𝒆 𝑧 + 𝒏 × 𝒗 𝑏 , where 𝒏 is the unitsurface normal vector directed into the fluid, 𝒗 is the fluid velocity, and 𝒗 𝑏 is the velocity of thesurface. Thus, the surface integrals in (2.63) and (2.64) can be re-written in terms of the vortexsheet strength and the body motion.We can easily develop discrete forms of the integrals (2.63) and (2.64) with the solutionsand notation described in this paper. For the volume integrals, let us denote diagonal matricescontaining the coordinates of the grid nodes by D x and D y . Thus, the expressions in (2.63) and(2.64) can be written in discrete form as 𝑃 𝑥 = Δ 𝑥 y 𝑇 w + r 𝑇𝑦 (cid:16) f + D n 𝑥 v 𝑏,𝑦 − D n 𝑦 v 𝑏,𝑥 (cid:17) , (2.66) 𝑃 𝑦 = − Δ 𝑥 x 𝑇 w − r 𝑇𝑥 (cid:16) f + D n 𝑥 v 𝑏,𝑦 − D n 𝑦 v 𝑏,𝑥 (cid:17) , (2.67)and Π = − Δ 𝑥 (cid:16) x 𝑇 D x + y 𝑇 D y (cid:17) w − (cid:16) r 𝑇𝑥 D r 𝑥 + r 𝑇𝑦 D r 𝑦 (cid:17) (cid:16) f + D n 𝑥 v 𝑏,𝑦 − D n 𝑦 v 𝑏,𝑥 (cid:17) . (2.68)8 − − − − − − − 𝑥 / 𝑅 𝑦 / 𝑅 − − − − − . . . 𝑋 𝑝 / 𝑅 𝑃 𝑥 𝑅 / Γ 𝑣 ( 𝑎 )( 𝑏 ) Figure 10. ( a ) Numerically simulated trajectories (——) of two point vortices ( • ) of opposite strengths Γ 𝑣 and − Γ 𝑣 being convected past a circular cylinder with radius 𝑅 , and ( b ) the 𝑥 component of the associated,numerically simulated impulse (——) in the fluid. Overlaid are the exact continuous trajectories (- - - -) andimpulse (- - - -). The simulation is performed with Δ 𝑥 / 𝑐 = . 𝛿𝑆 / Δ 𝑥 =
2, and Δ 𝑡 Γ 𝑣 / 𝑅 = . The overall force and moment exerted on the body are obtained from calculating these impulsesand computing their rates of change in (2.65). Part of this force and moment is attributable tothe dynamics of vorticity in the fluid. The remaining part is due to surface motion relative to thefluid, and we will discuss this in the next section.To illustrate the accuracy of the impulse-based calculation of force, we apply the method totwo examples. In the first example, we simulate the trajectories of two point vortices of oppositestrength, convected by each other past a cylinder. Figure 10 shows the trajectories and the 𝑥 component of the impulse together with their exact solutions. It is clear from the figure thatthe simulation shows very good agreement with the exact solution. In the second example, wecompare our simulation of the first instants of the unsteady, fully separated flow around a flatplate after impulsively starting a uniform flow with the Biot-Savart method from Darakananda &Eldredge (2019), using the same positioning rules to insert point vortices and the same time step.The vortex positions and the corresponding impulse and lift are compared in figure 11 and showvery good agreement as well.2.6.2. Added mass
The added mass tensor provides a measure of the inertial influence of the fluid on the bodyin response to changes in the body’s translational or rotational motion. The coefficients of theadded mass tensor of a body are obtained by computing the impulse components associated with9 − . . . − − . . 𝑥 / 𝑐 𝑦 / 𝑐 − . − . − . − − . 𝑃 𝑦 𝑐 / 𝑈 ∞ . . . . . . . . 𝑡𝑈 ∞ / 𝑐 𝐶 𝐿 ( 𝑎 ) ( 𝑏 )( 𝑐 ) Figure 11. Comparison of the simulated vortex shedding behind a flat plate of length 𝑐 at 60 ° in a uniformflow using the method in this paper ( ◦ and ——) and using the Biot-Savart method of Darakananda &Eldredge (2019) ( (cid:3) and - - - -). ( a ) The positions of the shedded point vortices, ( b ) the vertical componentof the associated impulse, and ( c ) the lift coefficient, one convective time after impulsively starting theuniform flow. At both edges, the Kutta condition is enforced. The simulation is performed with Δ 𝑥 / 𝑐 = . 𝛿𝑆 / Δ 𝑥 =
2, and Δ 𝑡𝑈 ∞ / 𝑐 = . a unit-valued component of motion (Eldredge 2019). The motion’s influence is both direct, viathe surface velocity, and indirect, in the bound vortex sheet that develops on the surface.For example, suppose that we consider translation at unit velocity in the 𝑥 direction, for whichthe motion is described by v 𝑏,𝑥 = , v 𝑏,𝑦 =
0, and s (cid:48) 𝑏 = r 𝑦 , and the associated bound vortexsheet—obtained without Kutta condition by solving the basic problem (2.12)—is f = S − P Γ s (cid:48) 𝑏 = S − P Γ r 𝑦 . The added mass coefficients corresponding to this motion are derived by substitutingthese into the impulse formulas (2.66)–(2.68): 𝑃 ( 𝑥 ) 𝑥 = r 𝑇𝑦 (cid:16) S − P Γ r 𝑦 − D n 𝑦 (cid:17) (2.69) 𝑃 ( 𝑥 ) 𝑦 = − r 𝑇𝑥 (cid:16) S − P Γ r 𝑦 − D n 𝑦 (cid:17) , (2.70) Π ( 𝑥 ) = − (cid:16) r 𝑇𝑥 D r 𝑥 + r 𝑇𝑦 D r 𝑦 (cid:17) (cid:16) S − P Γ r 𝑦 − D n 𝑦 (cid:17) . (2.71)Thus, the components of the added mass coefficients tensor associated with translation in the 𝑥 direction are 𝑚 𝐹𝑥𝑥 = 𝜌𝑃 ( 𝑥 ) 𝑥 (2.72) 𝑚 𝐹𝑥𝑦 = 𝜌𝑃 ( 𝑥 ) 𝑦 , (2.73) 𝑚 𝑀𝑥 = 𝜌 (cid:16) Π ( 𝑥 ) − 𝑿 𝑐 × 𝑷 ( 𝑥 ) (cid:17) , (2.74)where 𝑿 𝑐 is the centroid of the body, which can be calculated using (B 9), and the superscript 𝐹 and 𝑀 are used to denote the coefficient for the force and moment, respectively.A similar approach can be used to obtain the added mass coefficients due to unit translation inthe 𝑦 direction, for which v 𝑏,𝑥 = v 𝑏,𝑦 = , and s (cid:48) 𝑏 = − r 𝑥 . The coefficients due to unit rotationfollow from taking v 𝑏,𝑥 = − r 𝑦 , v 𝑏,𝑦 = r 𝑥 , and s (cid:48) 𝑏 = − ( D r 𝑥 r 𝑥 + D r 𝑦 r 𝑦 ) .0 − − − − 𝑥 / 𝑅 𝑦 / 𝑅 . . . . . . . 𝐺 / 𝑅 𝜇 / 𝛼 ( 𝑎 ) ( 𝑏 ) Figure 12. ( a ) Contours of the discrete streamfunction for an array of nine circular cylinders with radius 𝑅 ,spaced with a gap distance 𝐺 between each cylinder, of which the bottom left cylinder translates horizontally.( b ) Numerically simulated variation ( ◦ ) of the ratio between the principal value of added mass coefficienttensor and the largest self-added mass coefficient with the gap-to-radius ratio 𝐺 / 𝑅 . Overlaid are the values( (cid:3) ) obtained by Chen (1975) from solving a system of truncated analytical expressions. The simulations areperformed with Δ 𝑥 / 𝑅 = .
05 and 𝛿𝑆 / Δ 𝑥 = Multiple bodies
The previous sections provided the formulations for potential flow with the presence of abody. The extension of these expressions to multiple bodies is straightforward and consists ofallocating partitions of f to the different bodies. The surface streamfunction has to be partitionedaccordingly, with the body motion streamfunction s 𝑏 containing the values for the discrete surfacepoints from all the bodies and s allocating a uniform value to each body. The system (2.12) canthen be solved for the streamfunction field without modification.If we want to enforce an edge condition on one or multiple bodies, each body with edgeconditions has to add a Lagrange multiplier to the system. To accommodate the changes, weredefine 𝑠 as a vector, containing these Lagrange multipliers and redefine as a matrix with the ( 𝑖, 𝑗 ) th entry equal to one if 𝑖 is an index of f that belongs to the 𝑗 th body with edge constraintsand zero otherwise. With these new definitions, the saddle-point systems (2.29), (2.43), and (2.52)can be reused without further modification.Figure 12 demonstrates an example of a potential flow model with an array of nine circularcylinders and compares the ratio of the principal value of added mass coefficient tensor and thelargest self-added mass coefficient of the system with the results of Chen (1975).
3. Conclusion
We presented a two-dimensional, grid-based potential flow solver for use in a point vortex modelto simulate a steady or unsteady, unbounded flow around a body, with possible edge conditions.In the development of this solver, we made use of two algebraic techniques that allow us to mimicthe analytical treatment of potential flows around sharp-edged bodies. The first technique is tofollow the immersed boundary projection method and introduce a Lagrange multiplier for theno-penetration constraint in the streamfunction Poisson equation, thus forming a saddle-pointsystem. We can then identify the Lagrange multiplier as a discrete version of the continuousstrength distribution of a bound vortex sheet on the body, a singular vorticity distribution that iscommonly used to represent bodies in analytical potential flow. The second algebraic technique is1to decompose this discrete bound vortex sheet strength for sharp-edged bodies into a singular andsmooth part. We can then add constraints on the elements of the smooth part that are located atthe edges to make the overall discrete bound vortex sheet well-behaved. This allows us to enforcethe Kutta condition in a way that is similar to analytical treatments of the Kutta condition and tointroduce a way of enforcing generalized edge conditions. In a steady flow, one such constraintcan be enforced at a time by adding it to the saddle point system. The Lagrange multiplier isthen the—previously arbitrary—uniform value that can be added to the surface streamfunction.In an unsteady flow, where new point vortices can be introduced in the flow, we can enforceconstraints on multiple edges. This requires adding the edge constraints and the constraint of zerototal circulation to the saddle point system with their Lagrange multipliers being the uniformvalue of the surface streamfunction and the strengths of the new vortex elements. Furthermore,we can leverage the concept of the discrete bound vortex sheet strength to create expressions forthe impulse in the flow around bodies and for the added mass matrix for a system of arbitrarilyshaped bodies.We found that our method can replicate the results of a Biot-Savart method for the unsteadyflow around a flat plate with arbitrary precision. This motivates the goal of implementing a three-dimensional version of this solver since the concepts of the immersed boundary projection methodand the enforcement of edge conditions through the multiplicative decomposition of the vortexsheet strength generalize to three dimensions. Such a three-dimensional, grid-based solver wouldrely on the vector potential and a vector treatment of the vorticity field, and it could potentiallymake significant cost improvements over similar three-dimensional Biot-Savart methods. It is alsoimportant to note that our method is only specific in its use of constrained forces for the immersedboundary projection method. The idea behind our method is not restricted to the discretizationtools we used in our implementation.
Acknowledgements
The support for this work by the US Air Force Office of Scientific Research (FA9550-18-1-0440) with programme manager Gregg Abate is gratefully acknowledged.
Appendix A. Solution of general saddle-point systems
A general block system (with positive semi-definite matrix A ) can be decomposed as follows: (cid:20) A B 𝑇 B −C (cid:21) = (cid:20) A B S (cid:21) (cid:20) I A − B 𝑇 I (cid:21) , (A 1)where S ≡ −C − B A − B 𝑇 (A 2)is the Schur complement of the matrix system and I is the identity. By this decomposition, we candevelop an algorithm for the solution of the block system (cid:20) A B 𝑇 B −C (cid:21) (cid:18) xy (cid:19) = (cid:18) r r (cid:19) . (A 3)We will refer to x as the solution vector and y as the constraint force. We define the intermediatesolution vector ( x ∗ , y ∗ ) 𝑇 as the solution of the lower-triangular system (cid:20) A B S (cid:21) (cid:18) x ∗ y ∗ (cid:19) = (cid:18) r r (cid:19) (A 4)2and then the solution we seek can be found by back substitution of (cid:20) I A − B 𝑇 I (cid:21) (cid:18) xy (cid:19) = (cid:18) x ∗ y ∗ (cid:19) (A 5)The algorithm we derive from this is A x ∗ = r , S y ∗ = r − B x ∗ , (A 6) y = y ∗ , x = x ∗ − A − B 𝑇 y . It is also useful to have an inverse representation of the block matrix system: (cid:18) xy (cid:19) = (cid:20) A − + A − B 𝑇 S − B A − −A − B 𝑇 S − −S − B A − S − (cid:21) (cid:18) r r (cid:19) . (A 7) Appendix B. Some geometric relations for discrete surfaces
Consider a closed surface S 𝑏 with unit normal 𝒏 . We will recall some basic geometric relationshere, and then provide some discrete versions of these relations based on the set of points withcoordinates r 𝑥 , r 𝑦 , and and normal components n 𝑥 and n 𝑦 (which, the reader will recall, containthe surface length or area of each segment or panel associated with the points).The volume V 𝑏 of the region enclosed by S 𝑏 can be computed from the integral V 𝑏 = 𝑛 𝑑 ∫ S 𝑏 𝒙 · 𝒏 d 𝑆, (B 1)where 𝑛 𝑑 is the number of spatial dimensions (2 or 3). Using the notation above, the approximateform of this expression is V 𝑏 ≈ 𝑛 𝑑 ∑︁ 𝑗 r 𝑇𝑗 n 𝑗 , (B 2)where the sum is taken over the 𝑛 𝑑 components.An alternative formula for the volume is V 𝑏 𝒆 𝑗 = − 𝑛 𝑑 − ∫ S 𝑏 𝒙 × ( 𝒏 × 𝒆 𝑗 ) d 𝑆. (B 3)The components of this integral can be written discretely as V 𝑏 𝒆 𝑗 ≈ 𝒆 𝑗 𝑛 𝑑 − ∑︁ 𝑘 ≠ 𝑗 r 𝑇𝑘 n 𝑘 (B 4)And finally, a third alternative is V 𝑏 𝑰 = ∫ S 𝑏 𝒙𝒏 d 𝑆, (B 5)where 𝑰 is the identity. The discrete form of this is a diagonal matrix with r 𝑇𝑥 n 𝑥 , r 𝑇𝑦 n 𝑦 , and r 𝑇𝑧 n 𝑧 along the diagonal.Thus, we can conclude that the volume of the body is approximately V 𝑏 ≈ r 𝑇𝑥 n 𝑥 ≈ r 𝑇𝑦 n 𝑦 ≈ r 𝑇𝑧 n 𝑧 , (B 6)or any average of some combination of these.3The centroid of the body can be derived from the equation 𝑿 𝑐 V 𝑏 = ∫ S 𝑏 𝒙 · 𝒙𝒏 d 𝑆, (B 7) 𝑿 𝑐 × V 𝑏 𝒆 𝑗 = − 𝑛 𝑑 ∫ S 𝑏 𝒙 × (cid:2) 𝒙 × ( 𝒏 × 𝒆 𝑗 ) (cid:3) d 𝑆, (B 8) 𝑋 𝑐 ≈ V 𝑏 (cid:16) r 𝑇𝑥 D r 𝑥 + r 𝑇𝑦 D r 𝑦 (cid:17) n 𝑥 , 𝑌 𝑐 ≈ V 𝑏 (cid:16) r 𝑇𝑥 D r 𝑥 + r 𝑇𝑦 D r 𝑦 (cid:17) n 𝑦 (B 9) REFERENCESBaker, G. R. 1979 The “cloud in cell” technique applied to the roll up of vortex sheets.
Journal ofComputational Physics (1), 76–95.Beckers, D. & Eldredge, J. D. 2021 JuliaIBPM/GridPotentialFlow.jl v0.1.0. Available at: https://github.com/JuliaIBPM/GridPotentialFlow.jl/tree/v0.1.0 .Benzi, M., Golub, G. H. & Liesen, J. 2005 Numerical solution of saddle point problems. Acta Numerica , 1–137.Chatelain, P., Curioni, A., Bergdorf, M., Rossinelli, D., Andreoni, W. & Koumoutsakos, P. 2008Billion vortex particle direct numerical simulations of aircraft wakes. Computer Methods in AppliedMechanics and Engineering (13-16), 1296–1304.Chatelin, R. & Poncet, P. 2014 Hybrid grid-particle methods and Penalization: A Sherman-Morrison-Woodbury approach to compute 3D viscous flows using FFT.
Journal of Computational Physics ,314–328.Chen, S. S. 1975 Vibration of nuclear fuel bundles.
Nuclear Engineering and Design (3), 399–422.Chorin, A. J. & Bernard, P. S. 1973 Discretization of a vortex sheet, with an example of roll-up. Journalof Computational Physics (3), 423–429.Christiansen, J.P. 1973 Numerical simulation of hydrodynamics by the method of point vortices. Journalof Computational Physics (3), 363–379.Colonius, T. & Taira, K. 2008 A fast immersed boundary method using a nullspace approach and multi-domain far-field boundary conditions. Computer Methods in Applied Mechanics and Engineering (25-28), 2131–2146.Coquerelle, M. & Cottet, G.-H. 2008 A vortex level set method for the two-way coupling of anincompressible fluid with colliding rigid bodies.
Journal of Computational Physics (21), 9121–9137.Cottet, G.-H. & Koumoutsakos, P. 2000
Vortex Methods: Theory and Practice . Cambridge, UK:Cambridge University Press.Cottet, G.-H. & Poncet, P. 2004 Advances in direct numerical simulations of 3D wall-bounded flows byVortex-in-Cell methods.
Journal of Computational Physics (1), 136–158.Couët, B., Buneman, O. & Leonard, A. 1981 Simulation of three-dimensional incompressible flows witha vortex-in-cell method.
Journal of Computational Physics (2), 305–328.Cserti, J. 2000 Application of the lattice green’s function for calculating the resistance of an infinite networkof resistors. Am. J. Phys , 896–906.Darakananda, D. & Eldredge, J. D. 2019 A versatile taxonomy of low-dimensional vortex models forunsteady aerodynamics. J. Fluid Mech. , 917–948.Ebiana, A.B. & Bartholomew, R.W. 1996 Design considerations for numerical filters used in Vortex-in-cellalgorithms.
Computers & Fluids (1), 61–75.Eldredge, J. D. 2019 Mathematical Modeling of Unsteady Inviscid Flows , Interdisciplinary AppliedMathematics , vol. 50. Springer.Gazzola, M., Chatelain, P., van Rees, W. M. & Koumoutsakos, P. 2011 Simulations of single andmultiple swimmers with non-divergence free deforming geometries.
Journal of ComputationalPhysics (19), 7093–7114.Gillis, T., Winckelmans, G. & Chatelain, P. 2018 Fast immersed interface Poisson solver for 3Dunbounded problems around arbitrary geometries.
Journal of Computational Physics , 403–416. Hejlesen, M. M., Koumoutsakos, P., Leonard, A. & Walther, J. H. 2015 Iterative Brinkman penalizationfor remeshed vortex methods.
Journal of Computational Physics , 547–562.Katsura, S. & Inawashiro, S. 1971 Lattice Green’s functions for the rectangular and the square lattices atarbitrary points.
J. Math. Phys. , 1622–1630.LeVeque, R.J. & Li, Z. 1994 The Immersed Interface Method for Elliptic Equations with DiscontinuousCoefficients and Singular Sources. SIAM Journal on Numerical Analysis (4), 1019–1044.Liska, S. & Colonius, T. 2014 A parallel fast multipole method for elliptic difference equations. J. Comput.Phys. , 76–91.Marichal, Y., Chatelain, P. & Winckelmans, G. 2014 An immersed interface solver for the 2-Dunbounded Poisson equation and its application to potential flow.
Computers and Fluids , 76–86.Meng, J. C. S. & Thomson, J. A. L. 1978 Numerical studies of some nonlinear hydrodynamic problems bydiscrete vortex element methods. Journal of Fluid Mechanics (3), 433–453.Monaghan, J. J. 1985 Extrapolating B splines for interpolation. Journal of Computational Physics (2),253–262.Poncet, P. 2009 Analysis of an immersed boundary method for three-dimensional flows in vorticityformulation. Journal of Computational Physics (19), 7268–7288.Ramesh, K., Gopalarathnam, A., Granlund, K., OL, M. V. & Edwards, J. R. 2014 Discrete-vortexmethod with novel shedding criterion for unsteady aerofoil flows with intermittent leading-edgevortex shedding.
J. Fluid Mech. , 500–548.Rasmussen, J. T., Cottet, G.-H. & Walther, J. H. 2011 A multiresolution remeshed Vortex-In-Cellalgorithm using patches.
Journal of Computational Physics (17), 6742–6755.Rossinelli, D., Bergdorf, M., Cottet, G.-H. & Koumoutsakos, P. 2010 GPU accelerated simulations ofbluff body flows using vortex particle methods.
Journal of Computational Physics (9), 3316–3333.Spietz, H. J., Hejlesen, M. M. & Walther, J. H. 2017 Iterative Brinkman penalization for simulation ofimpulsively started flow past a sphere and a circular disc.
Journal of Computational Physics ,261–274.Wiegmann, Andreas & Bube, Kenneth P 2000 The Explicit-Jump Immersed Interface Method: FiniteDifference Methods for PDEs with Piecewise Smooth Solutions.
SIAM Journal on Numerical Analysis37