Accuracy comparison of different approaches for vortex sheet discretization on the airfoil in vortex particles method

Meshless MethodsAccuracy comparison of different approaches for vortex sheet discretization on the airfoil in vortex particles method V International Conference on Particle-based Methods – Fundamentals and Applications PARTICLES 2017 P. Wriggers, M. Bischoff, E. O˜nate, D.R.J. Owen, & T. Zohdi (Eds) ACCURACY COMPARISON OF DIFFERENT APPROACHES FOR VORTEX SHEET DISCRETIZATION ON THE AIRFOIL IN VORTEX PARTICLES METHOD KSENIIA S. KUZMINA1, ILIA K. MARCHEVSKY1, DARIO MILANI2AND EVGENIYA P. RYATINA1 1Bauman Moscow State Technical University (BMSTU), Applied mathematics department 2-nd Baumanskaya st., 5/1, 105005, Moscow, Russia E-mail: [email protected], iliamarc[email protected], [email protected] 2Bauhaus Universit¨at Weimar, Civil Engineering faculty Marienstrasse 13a, 99423 Weimar E-mail: [email protected] Key words: Vortex Particles Method, Boundary Integral Equation, Galerkin Approach, Finite Element Method, Discontinuity Extraction Abstract. The hierarchy of numerical schemes is developed for the solution of the boundary integral equation which arises in the vortex particles method. Two different approaches are considered and numerical schemes are presented, which provide 1-st and 2-nd order of accuracy with respect to two different norms. The developed finite element-type schemes have nearly the same computational cost as the simplest scheme with piecewise-constant solution. The necessary formulae are presented for linear system coefficients for the considered schemes. 1 Introduction One of the key questions which arise in implementations of meshless Lagrangian vortex particle methods for viscous incompressible flow simulation around the airfoil, is the choice of numerical approach for no-slip (or no-through in case of inviscid flow simulation) boundary condition satisfaction on the surface line of the airfoil. The airfoil in the flow in the most general case can be replaced with vortex sheets (attached and free) and attached source sheet. In coupled hydroelastic problems intensities of all these sheets are unknown; if law of motion for the airfoil is known, or it is rigid and immovable, only free vortex sheet intensity is unknown. The boundary condition according to one of two possible approaches (which will be described later) to its satisfaction can be written down in form of integral equation of the first or the second kind [1]. The kernels of such equations have very different properties. 1 691 Kseniia S. Kuzmina, Ilia K. Marchevsky, Dario Milani and Evgeniya P. Ryatina In the present paper the surface line of the airfoil is assumed to be approximated by polygon, which consists of rectilinear segments (usually called “panels”). The unknown vortex sheet intensity is approximated by piecewise-constant or piecewise-linear distribution. In the last case, it can be continuous or discontinuous, or continuous everywhere except some specified points that correspond to sharp edges and angle points of the airfoil. Some numerical schemes are considered and estimations for their accuracy and computational cost are obtained. 2 Brief description of vortex methods Vortex particles methods of flow simulation around airfoils are based on considering the vorticity as a primary computed variable and on the fundamental principle, which was discovered by prof. N.E. Zhukovsky in 1906: an immovable airfoil influences the inviscid incompressible flow just as attached vortex sheet placed on its surface line [2]. So it is possible to replace airfoil with vortex sheet and the main problem is how to determine the intensity of this vortex sheet. For some airfoil shapes (elliptical and Zhukovsky wing airfoils) it can be found by using conformal mappings technique [3], we use such solutions as benchmarks for verification of the numerical schemes being developed. When solving the Navier-Stokes equations the classical Zhukovsky approach remains correct in principle. However, it is required to account for the transfer of vorticity from the body surface to the flow which lead the vortex sheet to be free instead of being attached, Therefore, according to Lighthill’s approach [4] the phenomena can be modelled as a result of vorticity generation on the surface line Kdue to vorticity flux action within the small time period. Then the vorticity, which is concentrated in this vortex sheet with intensity γ(r), r∈K, becomes part of the vortex wake and moves in the flow according to the governing equations. When we consider a hydroelastic problems, — where either the law of motion is known a priori or coupled hydroelastic problem is being solved — in addition to free vortex sheet also attached vortex sheet (γatt) and attached source sheet (qatt) should be introduced [5, 6]. In the simplest case, when velocity of the surface line points VK(r) can be expressed explicitly as the function that depends on point position and time, the intensities of the attached sheets can be found explicitly: γatt(r)=VK(r)·τ(r),q att(r)=VK(r)·n(r),r∈K, where τ(r) and n(r) are tangent and normal unit vectors to the airfoil surface line, respectively, n(r)×τ(r)=k(kis unit vector orthogonal to the flow plane). Attached and free vortex sheet intensity vectors defined as γatt(ξ)=γatt(ξ)k,γ(ξ)=γ(ξ)k, are then introduced to reconstruct the velocity field in flow domain S, by employing the generalized Biot — Savart law [7]: 2 692 Kseniia S. Kuzmina, Ilia K. Marchevsky, Dario Milani and Evgeniya P. Ryatina V(r)=V∞+1 2πS Ω(ξ)×(r−ξ) |r−ξ|2dSξ+1 2πK γ(ξ)×(r−ξ) |r−ξ|2dlξ+ +1 2πK γatt(ξ)×(r−ξ) |r−ξ|2dlξ+1 2πK qatt(ξ)(r−ξ) |r−ξ|2dlξ,(1) where Ω(ξ)=Ω(ξ)k= curl Vis vorticity distribution in the flow domain. The crucial point is that the computation of the limit value of flow velocity on the surface line of the airfoil directly by using (1) is impossible. In fact, velocity field V(r) has jump discontinuity there [8]: vortex and source sheets from mathematical point of view are the curves at which tangent and normal components of velocity have jumps, respectively. In order to take them into account, the corresponding non-integral terms should be considered as elaboration of expression (1): V−(r)=V(r)−γ(r)−γatt(r) 2τ(r)+qatt(r) 2n(r),r∈K. (2) 3 Integral equation for vortex sheet intensity distribution Attached vortex sheet intensity distribution can be found from boundary integral equation which follows from boundary condition satisfaction on the surface line: V−(r)=VK(r),r∈K. (3) It can be shown [1], that there are two equivalent sufficient conditions for equality (3) satisfaction: •equality between normal components of velocities (N-scheme): V−(r)·n(r)=VK(r)·n(r),r∈K, (4) •equality between tangent components of velocities (T-scheme): V−(r)·τ(r)=VK(r)·τ(r),r∈K. (5) It is required to notice at this point that the resulting integral equations with respect to unknown γ(r), which can be derived from (4) and (5) after substituting there (2), have quite different properties. Usually in vortex methods boundary condition is being written down in form (4) and it leads to the 1-st kind singular integral equation: 1 2πK Qn(r,ξ)γ(ξ)dlξ=fn(r),r∈K, (6) where the kernel Qn(r,ξ) is unbounded (singular) of Hilbert-type Qn(r,ξ)=k×(r−ξ) |r−ξ|2·n(r)=−τ(r)·(r−ξ) |r−ξ|2. 3 693 Kseniia S. Kuzmina, Ilia K. Marchevsky, Dario Milani and Evgeniya P. Ryatina The right-hand side fn(r) depends on the airfoil shape, incident flow velocity, airfoil surface line velocity and vorticity distribution in the flow domain: fn(r)=−1 2qatt(r)−n(r)·V∞−VK(r)+ 1 2πS Ω(ξ)×(r−ξ) |r−ξ|2dSξ+ +1 2πK γatt(ξ)×(r−ξ) |r−ξ|2dlξ+1 2πK qatt(ξ)(r−ξ) |r−ξ|2dlξ.(7) For approximate numerical solution of equation (6) specific quadrature formulae are required, it can be discrete vortices-type quadrature formulae [8], which permits to pick out the principal value in Cauchy sense of the corresponding singular integral. When the numerical scheme of vortex method is being constructed on the basis of “tangent” approach (5), the resulting integral equation is the 2-nd kind Fredholm-type: 1 2πK Qτ(r,ξ)γ(ξ)dlξ−1 2γ(r)=fτ(r),r∈K, (8) where the kernel Qτ(r,ξ)=k×(r−ξ) |r−ξ|2·τ(r)=n(r)·(r−ξ) |r−ξ|2 is uniformly bounded function when Kis smooth curve and unbounded only in proximity to angle points or sharp edges of non-smooth airfoil surface line K, and right-hand side has the following form: fτ(r)=−1 2γatt(r)−τ(r)·V∞−VK(r)+ 1 2πS Ω(ξ)×(r−ξ) |r−ξ|2dSξ+ +1 2πK γatt(ξ)×(r−ξ) |r−ξ|2dlξ+1 2πK qatt(ξ)(r−ξ) |r−ξ|2dlξ.(9) It should be noted, that both equations (6) and (8) have infinite set of solutions. In order to select the unique solution, the additional equation should be added, which usually is being written down in the following form: K γ(ξ)dlξ=Γ.(10) The assignment of the value of the integral from solution over the surface line is the most common type for the additional condition for unsteady problems. 4 694 Kseniia S. Kuzmina, Ilia K. Marchevsky, Dario Milani and Evgeniya P. Ryatina 4 Problem discretization When solving problem using vortex methods, vortex wake (vorticity distribution in the flow) is normally being simulated by discrete vortex-type singularities — vortex elements: Ω(r,t)= n  w=1 Γw˜ δ(r−rw).(11) Here nis number or vortex elements; Γwand rware intensities and positions of vortex elements, respectively; ˜ δis two-dimensional Dirac delta function. Taking into account (11), we obtain for the integral term, which corresponds to the integral over the flow domain in right-hand side functions (7) and (9) 1 2πS Ω(ξ)×(r−ξ) |r−ξ|2dSξ= n  w=1 Γw 2π k×(r−rw) |r−rw|2. The simplest way to discretize the surface line Kof the airfoil is to approximate it by polygon consisting of Nrectilinear legs Ki, which we call hereinafter “panels”. Let’s denote lengthes of the panels as Li, their normal and tangent unit vectors as niand τi, respectively, i=1, ..., N. Now all the integrals over the surface line in the previous formulae can be written down as sums of integrals over the panels. 5 Numerical schemes for vortex sheet intensity computation For numerical solution of the integral equations (6) or (8) with additional condition (10) the ideas of the Galerkin method are used. According to this approach vortex sheet intensity distribution can be expressed as linear combination of some basis functions φk(r) with unknown coefficients, which can be found from the orthogonality condition of the residual p(r) of the corresponding integral equation to projection functions ψk(r):  K p(r)ψk(r)dlr=0.(12) Number of projection functions should be equal to number of unknown coefficients. For simplicity in this paper we suppose that there is no vorticity in the flow and the airfoil is immovable and non-deformable, so the intensities of the attached vortex and source sheets are equal to zero. These assumptions do not however limit the generality of the following considerations. The discretization procedure of the corresponding terms of the right-hand side (7) or (9) should be constructed similarly to free vortex sheet discretization. In this section we describe several numerical schemes, which follow from the different ways of the basis and projection functions choice. Four of them are based on T-scheme and one scheme corresponds to N-scheme. In order to estimate the accuracy of the numerical schemes the measures used are 5 695 Kseniia S. Kuzmina, Ilia K. Marchevsky, Dario Milani and Evgeniya P. Ryatina •a measure ∆γCh, which allows the average values ¯γiof computed vortex sheet intensity over the panels to be compared with average values of exact solution ¯γ∗ i: ∆γCh= max i|¯γi−¯γ∗ i|; •the traditional L1norm of the error: ∆γL1= K |γ(s)−γ∗(s)|J(s)ds ≈ N  k=1 sk+1  sk |γ(s)−γ∗(s)|Jkds, where γ(s) is approximate solution, parameterized by arbitrary variable swhose change corresponds to motion along the panels; J(s) is the Jakobian of such parametrization and Jkis its average value over the k-th panel; γ∗(s) means the projection of the exact solution to the corresponding point of the panel; sjis parameter value which corresponds to the beginning of the j-th panel. In all schemes the resulting system of linear algebraic equations is dense and nonsymmetric, so the most suitable way to solve it is Gaussian elimination method. 5.1 Scheme with piecewise-constant solution Introducing piecewise-constant basis functions φi 0(r), i=1, ..., N: φi 0(r)=1,r∈Ki, 0,r/∈Ki, we construct solution as piecewise-constant function γ(r)= N  i=1 γiφi 0(r). In this case the residual of integral equation (8) has the following form: p(r)= 1 2π N  j=1 γj Kj Qτ(r,ξ)dlξ−1 2 N  i=1 γiφi 0(r)−fτ(r). Choosing projection functions equal to the basis functions: ψi(r)=φi 0(r), we obtain the following system of linear algebraic equations from condition (12): 1 2π KiN  j=1 γj Kj Qτ(r,ξ)dlξdlr−1 2γiLi= Ki fτ(r)dlr,i=1, ... ,N, (13) with additional equation which follows from approximation of (10): N  i=1 γiLi=Γ.(14) 6 696 Kseniia S. Kuzmina, Ilia K. Marchevsky, Dario Milani and Evgeniya P. Ryatina The linear system (13), (14) is overdetermined; in order to regularize it the approach suggested by I.K. Lifanov [8] is used. The resulting system has the following matrix form: AI N L0γ R=b Γ,(15) where Ais square matrix block of size N×N; Aij =1 2π Ki Kj Qτ(r,ξ)dlξdlr−Li 2˜ δij, i, j =1, ..., N; ˜ δij is Kronecker delta; INis column consists of ones; L=(L1, ..., L N) is row consists of panel lengthes; bis right-hand side vector with components bi= Ki fτ(r)dlr,i=1, ..., N. γ=(γ1, ..., γ N)Tis vector of unknown coefficients; Ris regularization variable. Matrix of this system has size (N+ 1) ×(N+ 1), so Gaussian elimination procedure has computational cost (N+ 1)3/3≈N3/3. Piecewise-constant approximation provides the 2-nd order of accuracy for average values of vortex sheet intensity over the panels and only the 1-st order of accuracy in L1norm: ∆γCh∼N−2,∆γL1∼N−1. 5.2 Scheme with discontinuous piecewise-linear solution The simplest way to improve accuracy is by piecewise-linear solution representation over the panels as to substitute the piecewise-constant one. In addition to the earlier mentioned piecewise-constant basis functions φi 0(r), piecewiselinear basis functions φi 1(r) can be introduced, that vary from the value −1 2to 1 2over the i-th panel and are equal to 0 over the other panels: φi 1(r)=   (r−ci)·τi Li ,r∈Ki, 0,r/∈Ki, where ciis center of the i-th panel. Thus, the numerical solution is piecewise-linear discontinuous function which has the following form: γ(r)= N  i=1γiφi 0(r)+δiφi 1(r), where γiand δiare unknown coefficients. Residual of the integral equation (8) p(r)= 1 2π N  j=1 γj Kj Qτ(r,ξ)dlξ+δj Kj Qτ(r,ξ)φj 1(ξ)dlξ−1 2γ(r)−fτ(r) 7 697 Kseniia S. Kuzmina, Ilia K. Marchevsky, Dario Milani and Evgeniya P. Ryatina should be orthogonal to projection functions, which again are chosen equal to basis functions: ψk(r)=φi 0(r),φ i 1(r)N i=1. Total vorticity over the i-th panel doesn’t depend on δicoefficient, because the average value of φi 1basis function is equal to zero, so the approximation of the additional condition (10) has the same form (14). After regularization procedure, which is carried out similarly to the previous case, the resulting linear system has size (2N+ 1) ×(2N+ 1) and can be written down in form   A00 A01 IN A10 A11 ON LO N0   γ δ R =  b0 b1 Γ ,(16) where Apq,p, q =0,1 are square matrix blocks of size N×N Apq ij =1 2π Ki Kj Qτ(r,ξ)φj q(ξ)dlξφi p(r)dlr−DpLi˜ δij ˜ δpq, p, q =0,1, i, j =1, ..., N; ˜ δij and ˜ δpq are Kronecker deltas; D0=−1 2,D1=−1 24 ;γ=(γ1, ..., γ N)Tand δ=(δ1, ..., δ N)Tare vectors of unknown coefficients; INis column consists of ones; ONis row/column consists of zeros; L=(L1, ..., L N) is row consists of panel lengthes; b0and b1are vectors which form together the right-hand side: bp i= Ki fτ(r)φi p(r)dlr,i=1, ..., N. This numerical scheme provides the 2-nd order of accuracy both for ∆γChand ∆γL1errors, but the computational cost of this scheme is (2N+ 1)3/3≈8N3/3 and it is 8 times higher than the previously considered piecewise-constant scheme. 5.3 Finite element type numerical scheme (tangent) Note that in the framework of the piecewise-linear scheme numerical solution was considered to be discontinuous at the every node of the surface line discretization. However, for smooth airfoils the exact solution for vortex sheet intensity is continuous. This fact makes it possible to implement principles of the finite element method (FEM) and consider linear shape functions ˆ φiof the 1-st order (fig. 1) as basis and projection functions (the governing equation is being solved in the framework of tangent approach (8)). Functions ˆ φican be expressed through the previously introduced functions φi 0and φi 1: ˆ φi(r)=     1 2φi−1 0(r)+φi−1 1(r),r∈Ki−1, 1 2φi 0(r)−φi 1(r),r∈Ki. (17) 8 698 Kseniia S. Kuzmina, Ilia K. Marchevsky, Dario Milani and Evgeniya P. Ryatina 1−i K i K 1+i K 1i− c i c 1i+ c 1() i− γr () i γr 1() i+ γr 1 ˆ() i φr 1−i K i K 1 ˆ() i+ φr 1 ˆ() i− φr Figure 1: Linear shape functions The resulting regularized linear algebraic system has the form ˆ AI N ˆ L0ˆγ R=ˆ b Γ,(18) where ˆγ= (ˆγ1, ..., ˆγN)Tis vector of unknown variables; ˆγiis the value of vortex sheet intensity at the i-th node (which coincides with the beginning of the i-th panel). Coefficients of matrix block ˆ A(which size is N×N) can be expressed through the previously introduced coefficients of matrices Apq,(p, q =0,1): ˆ Aij =1 4A00 ij +A00 i−1,j +A00 i,j−1+A00 i−1,j−1+1 2−A01 ij −A01 i−1,j +A01 i,j−1+A01 i−1,j−1+ +1 2−A10 ij +A10 i−1,j −A10 i,j−1+A10 i−1,j−1+A11 ij −A11 i−1,j −A11 i,j−1+A11 i−1,j−1. Components of vector ˆ Land right-hand side vector ˆ balso can be expressed through the coefficients introduced in section 5.2: ˆ Li=1 2Li+Li−1,ˆ bi=1 2b0 i+b0 i−1−b1 i−b1 i−1. The matrix size now is (N+1)×(N+ 1), so the computational cost of the solution procedure remains the same as for piecewise-constant scheme (N3/3), but it provides the 2-nd order of accuracy for both errors ∆γChand ∆γL1in case of flow simulation around smooth airfoils. If the surface line of the airfoil has sharp edges or angle points, the exact solution is discontinuous and can be unbounded in proximity to these points, so the developed approach can’t be applied. 5.4 Finite element type numerical scheme (normal) This section presents a common approach widely employed in bluff body aerodynamics e.g. in bridge aerodynamics [9], which is the counterpart of the scheme devised in section 5.3. The present scheme considers the “normal” approach framework to enforce the boundary conditions. In order to solve numerically the integral equation (6), basis functions are chosen the same (linear shape functions (17)), but the projection functions now are Dirac delta functions: ψi(r)=˜ δ(r−ci), where ciis center of the i-th panel. 9 699