CC BY
18Mathematics i Mathematical Modelling
Mathematics and Mathematical Modeling, 2018, no. 6, pp. 22-51.
DOI: 10.24108/mathm.0618.0000158
Received:
19.11.2018
Electronic journal http://mathmelpub.ru ISSN 2412-5911
© NEICON
Improved Algorithm of Boundary Integral Equation Approximation in 2D Vortex Method for Flow Simulation Around Curvilinear Airfoil
Kuzmina K. S.1'2, Marchevsky I. K.1'2 Soldatova I. A.1
1Bauman Moscow State Technical University, Moscow, Russia 2Ivannikov Institute for System Programming of the RAS, Moscow, Russia
The problem of the accuracy improving is considered for vortex method. The general Galerkin-type approach is considered for the numerical solution of the boundary integral equation. It is shown that the airfoil surface line representation as a polygon consisting of rectilinear panels can lead to incorrect behavior of the numerical solution of the boundary integral equation with respect to unknown vortex sheet intensity, especially in case of considerably different lengths of the neighboring panels. However, in this case, the exact analytical formulae are derived for the coefficients of the linear algebraic system, which approximate the integral equation. The approach is developed which makes it possible to take into account explicitly the curvature of the airfoil surface line. The approximate analytical formulae are derived for the coefficients of the linear system, which are expressed through the integrals of the fractions, which denominators can be equal to zero. Numerical example is considered for the elliptic airfoil in the incident flow at its essentially non-uniform discretization. The suggested algorithm has acceptable numerical complexity, however it permits to obtain rather good numerical results.
Keywords: Vortex method, boundary integral equation, 2D flow simulation, curvilinear panel, numerical integration, non-uniform discretization
Despite the fact that vortex methods [1,2] are known for more than 50 years, it is hardly possible to consider their theory to be finished [3, 4, 5]. Even for 2D implementations of vortex methods, which are much easier in comparison with 3D case, some questions are still urgent, both in "mathematical background" (i.e., proof of the solution existence and uniqueness, numerical solution convergence to the exact one, convergence rate estimation), and related to the methodology of the numerical schemes optimal choice for numerical simulation in different problems, their accuracy and numerical complexity estimations, at least on empirical or semi-empirical level, etc.
Hereinafter we consider only 2D problems of flow simulation for viscid (or with some "minor" modifications of the algorithm — also for inviscid) incompressible flow. Such problems are rather simple due to orthogonality of the vorticity vector Q = curl V to the flow plane, in particular, the scalar product (V • Q) is always equal to zero.
Introduction
1. Basic idea of the vortex methods
The main idea of the vortex methods is to solve the system consisting of the continuity equation and the Navier — Stokes equation
V- V = 0, (1)
+ (V ■ V)V = -Vp + vV2V (2)
ot p
with traditional notations for all the flow variables and physical properties of the medium by using pure Lagrangian (meshless) approach.
The boundary conditions are considered to be the following:
1) the perturbation decay conditions at the "unbounded" part of the computational domain:
V(r,t) ^ Vp(r,t) ^ p^, |r| > to, (3)
where p^ = const, and V^ is constant vector field;
2) the no-slip condition at the airfoil surface lines:
V(r,t) = Vk(r,t), r e K(t), (4)
where K (t) is the surface line of the airfoil, or set of the surface lines in case of system of airfoils being considered in the flow; VK(r, t) is the velocity of K(t) at the corresponding point.
We mention briefly the main points of the Viscous Vortex Domains (VVD) method developed by prof. G.Ya. Dynnikova and described in details in [6,7,8]. We introduce the particles which carry the vorticity, which, in turn, is considered as the primary computational variable. The incompressible flow velocity field can be reconstructed by using the generalized Helmholtz decomposition [9, 10, 11], which can be applied for both bounded or unbounded flows in case of the no-slip boundary conditions on the airfoil surface lines
V(r t) V + f Q(^,t) x (r - A dS + I Y(£,t) x (r - £) + V(r,t) = V- +/ 2n|r - ||2 dS* + / 2n|r - ||2 ^ +
S(t) K(t)
t (VK(i.i) ■ T(i,i))fc x (r - O J- (VK(i,i) ■ n(i,i))(r - i)
^ ^-^-^ ^—^—, (5)
K(t) 1 1 K (t) 1 1
where S(t) is the flow domain (or the part of the flow domain with non-zero vorticity, because the integration is performed over the vorticity); n is the outer with respect to the airfoil normal unit vector; t is the tangent unit vector (n(£, t) x t(£, t) = k; k is the unit vector orthogonal to the flow plane); y(£, t) at £ e K is the vortex sheet [12, 13] intensity at the airfoil surface line, which can be considered as the result of vorticity flux action during some time period, where vorticity flux is the vorticity which is generated at the unit length of the airfoil surface line in unit time.
Note that vorticity flux influences directly (algebraically) the time derivative of the flow velocity
dV
—7 so it seems to be not very suitable for practical computations, however it can be useful, firstly, for analytical description of the background of vortex methods, and secondly, for development of high-order numerical schemes in future perspective, not in the framework of the present work.
Tangent and normal components of the airfoil surface line velocity can be considered as intensity of the attached vortex and source sheets, respectively,
Vk(£,t) ■ t(£,t) = Yatt(l,t), Vk(£,t) ■ n(£,t) = qatt(£,t), £ e K,
so the formula (5) can be written down in the following way:
v (r,t)=v „+/n(tt;,x(£|-£) ds5+
S(t)
r (y(£,t) + Yatt(£,t)) X (r - £) r qatt(£, t)(r - £) J7 ^
+ f ---^ +f 2(;,^£|2£) ^, (6)
k (t) 1 K(t)
where y(£, t) = y(£, t)k, Yatt(£, t) = Yatt(£, t)k.
The formulae (5) and (6) remain correct for arbitrary point r in the flow domain, but at its boundary (i.e., at the airfoil surface line K), according to properties of the generalized Helmholtz decomposition [9, 11], the following equation is satisfied:
a(r, t) Vk(r, t) = VM + J £ - £) dS? - a(r, t)7(r, t)T(r, t) +
s(t) n|r £|
r (y(£,t) + Yatt(£,t)) X (r - £) r qatt(£,t)(r - £)
+ f '--+ $ 2n|r)--£|2£) , r e K, (7)
k (t) Sl K(t)
where a(r, t) = 1 for smooth parts of the airfoil surface line; for other cases a is equal to the outer angle between the tangent lines at the corresponding point, divided by 2n. The contour integrals in the right-hand side in (7) are singular, they should be considered in sense of the principal values [3, 14].
We consider the law of motion of the airfoil to be known, so the region S(t), the shape and the position of the surface lines K (t), and intensities of the attached sheets yatt(r, t) and qatt(r, t), r e K(t), are also known. This assumptions mean that either "classical" aero- or hydrodynamics problem is considered (flow simulation around immovable airfoil or around the airfoil, which moves according to the given low), or coupled fluid-structure interaction (FSI) problem in the framework of the partitioned approach with some coupling strategy. However, it can be generalized for the FSI problems being considered in the framework of monolithic approach [15], according to
which the law of motion of the airfoil, and consequently yatt(r,t) and qatt(r,t), r e K(t), are unknown. They can be expressed through the surface line velocity, which, in turn, is considered to be unknown together with vortex sheet intensity 7(r, t). The coupling scheme in this case is based on linear relationship between hydrodynamic force acting the airfoil, 7(r, t), Yatt(r, t) and qatt(r, t) [15, 16].
It is necessary to satisfy the equation (7), which follows from the no-slip boundary condition at the airfoil surface line. We consider it for the above described simplest case, when there is only one unknown 7(r, t). It can be written down in the form of vector boundary integral equation of the second kind:
/ Y(£,t) dk - a(r,t)Y(r,t)T(r,t) = f (r,t), (8)
where f(r, t) is known function, which depends on the incident flow velocity, the vorticity distribution in the flow domain and the law of motion of the airfoil surface line:
f (r, t) = a(r, t)VK(r, t) - VM - J ^Jr 2 £) dS? -
S(t) n r
2n|r - 5|2
Y-'ii.t) x <r - €) dk _ r q-^Xr - & r e K(t). (9)
2n|r - £|2 s J 2n|r - £
K(t) K(t)
The equation (8) has infinite set of solutions, in order to select the unique one, the additional condition should be added. Normally in vortex methods this condition is isoperimetric-type condition
f Y(r,t) dlr = r, (10)
K(t)
where r is the given value of total circulation around the airfoil. The condition of such type normally follows from "physical" problem statement.
Coming back to the governing equation, one can notice that now the incompressibility equation (1) is satisfied, as well as the perturbation decay condition (3) for velocity field. The no-slip condition (4) is satisfied through vortex sheet intensity choice according to (8)-(10). So, the last equation, which should be satisfied, is the Navier — Stokes equation (2), which for 2D flows can be written down in Helmholtz-type form (with respect to unknown Q = curl V instead of V and p):
dt- + Vx (ft x (V + W)) = 0, (11)
where W is the so-called diffusive velocity field [17], which takes into account the viscosity influence.
The equation (11) can be considered as the transfer equation for vorticity Q = Qfe, which moves along the streamlines of the vector field V + W. The last proposition makes it easy to solve
this problem by using Lagrangian meshless method: inside the flow domain we introduce particles (vortex elements, called in [7, 8, 15] "viscous vortex domains") at points r», i = 1, ..., n, which carry the vorticity. In order to solve the equation (11), it is enough to consider that the strengthes of the vortex elements r», i = 1, ..., n (which have sense of the velocity circulations around the vortex elements if they are considered as point vortices) remain constant, while their motion is described by the system of ordinary differential equations
^ = V + W, i = 1, ..., n, dt
which is solved, as a rule, by using the explicit Euler method, however, higher-order explicit methods also can be used [1].
If the vorticity distribution in the flow domain is known, the velocity field V can be reconstructed by using the above discussed Biot — Savart law (5), the efficient approach to diffusive velocity W computation is developed in [6, 18].
In order to reconstruct the pressure field, the analogue of the Cauchy — Lagrange integral can be used, which is suitable for viscous incompressible non-potential flows [18]. It follows from the equation (2), and the resulting pressure distribution satisfies the perturbation decay condition (3) for pressure. Moreover, it is possible to derive rather simple expressions for integral hydrodynamic force and torque, acting the airfoil in the flow [15].
So, the described pure Lagrangian meshless approach makes it possible to satisfy in 2D case the governing equations (1)-(2) and the boundary equations (3)-(4).
Numerical experiments show that in order to provide high accuracy, the boundary integral equation (8) solution becomes the "critical step".
2. The brief review of the known approaches to the boundary integral equation solving in vortex methods
There are two possible approaches to solve the equation (8), which is vector equation. It is mentioned in [3, 1] and shown [11] that it is possible to consider either scalar integral equation which corresponds to equality of normal components of flow velocity and the surface line velocity, or scalar equation which corresponds to equality of tangent velocity components. These approaches are equivalent from a mathematical point of view, however they lead to quite different equations.
2.1. The N-scheme
Multiplying the left and right-hand sides of the equation (8) by the unit normal vector n(r, t), we obtain an integral equation of the first kind
/ Qn(r,£)y(£,t) d/? = fn(r,t), r e K(t), (12)
K(t)
where
«Wr £) = k2X|r(r-"£]l) ■ = -T£), A«^ = f ■ "M)
and it is assumed, as earlier, that n(r, t) is the unit outer normal vector, and t(r, t) is the tangent unit vector, which direction is chosen such as n(r, t) x t(r, t) = k.
The equation (12) is the singular integral equation with Hilbert-type singular kernel Qn, the integral should be considered in Cauchy sense. Such equations are very common in the framework of boundary integral equations methods since they follows straightforwardly from the Newmann problem for the Laplacian equation, when the solution is represented as a double layer potential. The main difficulty in such problems is its approximation with high accuracy. Numerical experiments show that widely-used numerical schemes provide rather high error, especially in case when the right-hand side fn(r,t) has high gradients, which takes place, in particular, at flow simulations in unsteady problems, when flow separation is observed and there is non-zero vorticity in the boundary layer, which is represented by a set of vortex elements, placed in the neighborhood of the airfoil surface line [19].
2.2. The T-scheme
According to the alternative approach, suggested by S. Kempka and co-authors [11], one can multiply the left and right-hand sides of the equation (8) by unit tangent vector t(r,t). The resulting equation is the Fredholm-type integral equation of the second kind
/ qt(r,£)y(£, t) - a(r,t)Y(r,t) = f (r,t), r e K(t), (13)
k (t)
where
qt £) = ■ T <-> = , fT m) = f ■ T M)
It is important to note that the kernel QT (r, £) of the equation (13) is bounded for smooth airfoils; if the surface line is C2-smooth curve, it is easy to prove that
H*
lim 0 |qt(r,£)| = J-, r e K(t), £ e K(t),
where k* is the curvature of the curve K (t) at the corresponding point.
If there is angle point at the airfoil surface line K (t), than the kernel QT becomes unbounded and the solution has weak singularity in proximity to this point. For some particular cases the asymptotic behavior of the solution is described in [3].
The numerical schemes, based on the T-scheme, provide much higher accuracy in comparison with N-schemes, especially in the above mentioned case of flows simulation with vortex elements in the neighborhood of the airfoil surface line [20, 21, 22].
2.3. Numerical method of the boundary integral equation solution
Hereinafter we consider the particular case when there is one airfoil in the flow, however, all the results can be easily generalized to the case of flow simulation around multiple airfoils.
For the numerical solution of the boundary integral equation, either (12) or (13), the so-called panel methods [23] are used, when the initial airfoil surface line is replaced by panels (or boundary elements) and the integral equation, being approximated on these panels, is reduced to a system of algebraic equations. The most general way to approximate the boundary integral equation (12) or (13), is based on the application of Galerkin approach [24]. Let us consider the integral equation in the following form
/ Q(r, dlf - ^a(r,t)Y(r,t) = f(r,t), r e K(t), (14)
K(t)
where Q = Qn, f = fn, ^ = 0 for the equation (12) and Q = QT, f = fT, ^ =1 for the equation (13).
We assume now that the airfoil surface line of the airfoil is parameterized with the arc length, than the equation (14) has the form
L
J Q(s, o)y(a) do - pa(r,t)Y(s) = f (s), s e [0, L], (15)
0
where L is the length of the airfoil surface line; hereinafter we omit dependencies of the kernel, the right-hand side and the solution of the integral equation on time; such dependencies are not essential from the point of view of its numerical solution.
The following way to solve the equation (15) numerically seems to be the most efficient.
1. The surface line of the airfoil is split into N parts, usually called "panels", which endings correspond to the arc length parameter values Sj, i = 0, ..., N, where s0 = 0, sN = L; the i-th panel corresponds to the parameter values in range s e [si-i, sj].
2. The set of basis functions {^,q (s)}, i = 1, ..., N, q = 0, ..., m, is introduced; we assume that for the given i the functions ^i,q(s) can differ from zero only at the i-th panel.
3. The approximate solution of the equation (15) is expressed as the combination of the basis functions with unknown coefficients:
N m
y(s) = EE Yq ^q (s). (16)
i=1q=0
4. In order to find the values of the coefficients y] we require residual of the equation (15) after substitution (16) to be orthogonal to the set of projection functions {^(s)}, i = 1, ..., N,
p = 0, ..., m. We again consider ^f(s) = 0 at s e [si-1, si]:
N m Si Sj m
EEy] ^p(s)ds / Q(s,°M(o)do-Mr,t)EYq / (s)^p(s)ds =
j=1 q=0 Si-i s/_i q=0 Si-i
Si
J f (s)^P(s) ds, i = 1, ...,N, p = 0, ...,m. (17)
Si— 1
Note, that it is proved in [11], that in the framework of the T-scheme there is no need to specify any additional condition, i.e., the equation (13) has unique solution, as well as linear algebraic system (17) approximating it. However, this applies only for the case of the bounded flow domain; in the considered case of unbounded flow domain both integral equations (12) and (13) have infinite set of solutions, the unique solution can be selected by taking into account the additional condition (10), which under the considered assumptions has the following form:
N m
EEyÍ / tf(s) ds = r. (18)
i=l p=0 si—i
The resulting system (17)-(18) is overdetermined in general case since it consists of N • (m + 1) + 1 equations with respect to N • (m +1) unknown coefficients yp. It can be regularized by using the technique developed by prof. I.K. Lifanov [3].
Note that the implementation of the Galerkin approach with traditional shape functions, which refer not to the panels but to the nodes and have non-zero values at two neighboring panels, also can be expressed in terms of the above described algorithm. The description of the corresponding algorithm can be found in [20]. Known numerical schemes for numerical solution of the integral equation (8) in 2D vortex methods permit to consider unknown vortex sheet intensity distribution to be singular (basis functions are Dirac delta-functions) [3], piecewise-constant [21] and piecewise-linear [20, 22] along the panels on the surface line. Projection functions can be either the same as the basis ones, or different from them(that corresponds to the Galerkin — Petrov method).
3. The airfoil surface line discretization with rectilinear panels
The easiest way to approximate the surface line of the airfoil is to replace it with the polygon which vertices coincide with the panels endings; in the above introduced notations — points r(s¿), i = 0, ..., N. Moreover, in this case the normal and tangent unit vectors remain constant along the panels, note them n and tj, respectively. Denoting
1 r(s) - r(a)
G(s, a) —
2n |r(s) - r(a)|2 '
the gradient of the Green's function — fundamental solution of the 2D Laplacian equation, taking into account that in (9) a = - for all the panels (except the panels endings; this however doesn't effect the result of integration along the panels in the framework of Galerkin approach), it is possible to write down the equation (17) as follows:
N m si Sj
EE Y* e • / # (s) ds f G(s,a)^q (a) da -
j=1 q=0 si—i sj— i
- f E Yq / tf(s)#(s) ds = n ■ J (- V«)#(s) ds -
q=0 Si-1 Si-1
s; f n(£,t) X fr(s) -^ * s; 7
- ni ■ J dsj 2 ■r(s)- £|2 ^P(s) dSî - Eei ■ J $"(s) ds y G(s,a)Yatt(a) da -
s
-1 S(t) 2n|r(s) I1 j=l s/_i
- E n • j ^f(s) ds J G(s,a)qatt(a) da, i = 1, ... , N, p = 0, ..., m, (19)
j=1 Si-1 Sj-1
where the vectors ni3 ei and parameter - are the following:
• for the N-scheme ni = ni; ei = —ti; - = 0;
• for the T-scheme ni = ti; ei = ni; - = 1.
We assume for consistency that the intensities of the attached sheets Yatt(s) and qatt(s) are also expressed through the basis functions:
N m N m
Y att(s) = EE (y att)PV?(s), qatt(s) = EE(?att)M(s). (20)
i=1p=0 i=1p=0
We consider the vorticity to be represented by a set of vortex elements with circulations rw, placed at the points pw, w = 1, ..., n:
N
fi(r)= E rwi(r — Pw),
w=1
where £(r — pw) is the 2D Dirac delta-function.
Then the linear equations (19) can be written down in the following form:
N m si Sj — m si
EEYqei• J (s)ds J G(s,a)^q(a)da — -EYq J (s)^f(s)ds =
j=1 l=0 Si-1 Sj-1 l=0 Si-1
si t -,- / \ N si
ni ■ f ("T^" - V«)^P(s) ds - ei ■ EE rw j Gw(s)^f(s) ds -
o. , w = 1 „. ,
N m
- EE( (Yatt)q ei + (qatt)q nil X
j=iq=0
si ° J
X J ^P(s) ds J G(s,a)^q (a) da, i =1,...,N, p = 0, ..., m. (21)
Si-1 Sj_1
Here, Gw(s) = Ge(r(s) — pwJ is smoothed Green's function gradient, introduced to avoid unbounded velocities induced by vortex elements, e.g.
1 r(s) - Pw
GJr(s) - Pw
2n max! lr(s) - Pw|2, e2}
For the N-scheme, the value of the parameter e as a rule has the same order as panel size, for the T-schemes it is small positive value that is required only for division by zero avoiding.
Different families of functions can be used as basis and projection functions, for example, constant and linear functions can be considered:
(r - Ci) ■ Ti
, , Li
0 r f Ki; I 0, r / Ki.
o, ^ ) 1 r f Ki; ^ J —t-, r f Ki;
Here, Ki, i = 1, ..., N, is the i-th rectilinear panel; ci is the center of the i-th rectilinear panel; Li is the length of the i-th curvilinear panel. Put simply, the basis function ^ (r) varies linearly with respect to the arc length from the value - - at the beginning of the i-th panel to the value - at its ending.
Note that for constant or linear basis and projection functions (or Dirac delta-functions) in arbitrary combinations, all the integrals in (21) can be calculated exactly by using formulae given in [25].
The described approximation scheme seems to be rather simple and suitable for practice. However, it works well only in the cases of close to uniform airfoils surface line discretization. This fact will be illustrated below for the simplest test case.
We consider the elliptic airfoil with 2:1 semiaxes ratio (a =1, b = 0.5, its parametric equations x(p) = a cosp, y(p) = b sinp). The airfoil is immovable (VK = 0, Yatt = 0, qatt = 0), vortex wake is absent (n = 0), the incident flow velocity has unit magnitude | V= 1 and it is directed at an angle P = n/6. Numerical computations are performed for the T-scheme.
The exact solution for vortex sheet intensity can be easily obtained by using the complex potential theory and conformal mappings technique [26]:
a + b
7.(p) = 2M sin(p - p) ■ + b^ |. (22)
3.1. Quasi-uniform airfoil discretization
Firstly, we discretize the elliptic airfoil into N = 20 panels by the following way: the endings
2 nj
of the panels correspond to the parameter p values pj = , j = 0, ..., N. The positions of these points are shown in Figure 1. Panels lengths are not the same, but the lengthes of the neighboring panels differ slightly. Such discretization we call "quasi-uniform".
The results of numerical solution of the linear system (21), (18) according to the above described algorithm are shown in Figure 2 in comparison with the exact solution (22): a — for the case when numerical solution is assumed to be piecewise-constant (m = 0); b — for the piecewise-linear numerical solution (m = 1).
It is seen that both numerical solutions are "good", they are close to the exact solution. So the developed algorithm in this case is quite suitable for practical purposes, especially taking into account that for all the integrals in (21) there are exact analytical expressions derived in [25].
Figure 1: Quasi-uniform ellipse discretization into 20 panels
(a) (6)
Figure 2: Piecewise-constant (a) and piecewise-linear (b) numerical solution for the rectilinear panels at quasi-uniform ellipse discretization; red curve corresponds to the exact solution; vertical lines correspond to the panels endings
3.2. Non-uniform airfoil discretization
Now we consider essentially non-uniform airfoil discretization. The upper arc of the ellipse, as earlier, is split into 10 panels, but the lower arc we split now into 50 panels (Figure 3). The maximal ratio of the neighboring panels lengths now is about 5.2.
v
^^ ~ 0.4 / 0.2 x.........
1-1.0 -0.5 \ -0.2 0.5 1.01
Figure 3: Non-uniform ellipse discretization into 50 panels
(a) (6)
Figure 4: Piecewise-constant (a) and piecewise-linear (b) numerical solution for the rectilinear panels at non-uniform ellipse discretization; the red curve corresponds to the exact solution; the vertical lines correspond to the panels endings
Application of the above described algorithm (with rectilinear panels representation) leads to the results, shown in Figure 4, again, in comparison with the exact solution (22): a — for the case when numerical solution is assumed to be piecewise-constant (m = 0); b — for the piecewise-linear numerical solution (m = 1).
Now the piecewise-constant solution (Figure 4, a) remains qualitatively correct, however, there is significant error for the panels that adjust to the large semiaxis of the ellipse (that are in the neighborhood to the border between the panels of significantly different length — the left and right quadrant points). The piecewise-linear numerical solution in these regions is incorrect.
The parts of these pictures, but at more obvious scale, are shown in Figure 5 and 6.
(a) (6)
Figure 5: Piecewise-constant (a) and piecewise-linear (b) numerical solution for the rectilinear panels in the neighborhood of the left quadrant point of the ellipse at non-uniform ellipse discretization; the red curve corresponds to the exact solution; the vertical lines correspond to the panels endings
(a) (5)
Figure 6: Piecewise-constant (a) and piecewise-linear (b) numerical solution for the rectilinear panels in the neighborhood of the right quadrant point of the ellipse at non-uniform ellipse discretization; the red curve corresponds to the exact solution; the vertical lines correspond to the panels endings(solution is continued to the right according its periodicity)
Now we briefly explain the source of such an error of the numerical solution. The problem is connected with "dummy" angle points which arise at the airfoil discretization into rectilinear panels (Figure 7).
Original surface line
Figure 7: Surface line of the initial airfoil(dashed line) and its discretization into rectilinear panels with "dummy" angle point formation
In case of rectilinear discretization of the surface line, the exact solution for the polygonal airfoil has singularities at every angle point, and the more carefully it is represented by using some numerical method, the more significant difference it would have with the exact solution for the initial smooth airfoil. However, when the numerical solution is assumed to be piecewise-constant along the panels, or even piecewise-linear, but the neighboring panels have nearly the same lengths (Figure 1), airfoil approximation with the polygon is suitable (Figure 2). For significantly different lengths of neighboring panels (Figure 3), the result of numerical solution can be far from exact solution, not only quantitatively, but even qualitatively (Figure 4-6).
So, the mentioned problem can be solved only by taking into account the curvilinearity of the airfoil.
4. The airfoil surface line discretization with curvilinear panels
Let us take into account the curvilinearity of the airfoil explicitly, by numerical calculation of all the integrals in the left and right-hand sides of (17), considering it as earlier in the framework of the T-scheme. Note that in the considered particular case the right-hand side is determined only by the incident flow:
f (s) = -V• n(s) = -V• t(s).
However, we present the approximate expressions for all the other terms arising in the right-hand side of (17), which correspond to the influence of the vortices modelling the vortex wake, attached vortex and source sheets (for movable airfoil) and the surface line velocity.
4.1. Direct integrals computation
For computations the Wolfram Mathematica computer algebra system was used. The results are shown in Figure 8 and 9 only for the regions in the neighborhood of the left and right quadrant points of the ellipse.
It is seen that the solution now is rather "good", the nun-uniformity of the discretization practically doesn't effect the numerical solution. However, this method hardly can be used in practice, because the numerical calculation of the integrals is too complicated problem.
Now we suggest the efficient approach that on the one hand, permits to take into account the curvilinearity of the airfoil correctly and deal with non-uniform airfoil discretization, and on the other hand, has acceptable numerical complexity.
(a) (6)
Figure 8: Piecewise-constant (a) and piecewise-linear (b) numerical solution for the curvilinear panels in the neighborhood of the left quadrant point of the ellipse at non-uniform ellipse discretization; the integrals along the curvilinear panels are calculated numerically by using the Wolfram Mathematica software
(à) (6)
Figure 9: Piecewise-constant (a) and piecewise-linear (b) numerical solution for the curvilinear panels in the neighborhood of the right quadrant point of the ellipse at non-uniform ellipse discretization; the integrals along the curvilinear panels are calculated numerically by using the Wolfram Mathematica software
4.2. Approximate left-hand side integrals computation
We consider piecewise-linear numerical solution (m = 1), so as for rectilinear panels, two families of basis functions are introduced:
s(r) - ) r e K.
H>\ (r ) = \ * ' "
0, r / K.
^0(r) =
1, r G Ki ; 0, r G Ki;
Here, Ki, i = 1, ..., N, is the i-th curvilinear panel; s(r) is the airfoil surface line arc length function (natural parameter) value at the point r; c is the center of the i-th curvilinear panel; Si is the length of the i-th curvilinear panel.
Projection functions are considered to be equal to basis functions, ^p = ^p, i = 1, ..., N,
p = 0, 1.
The approximate solution has the following form:
N
Y (r ) = E (Y0 ^0 (r )+ Y1 (r )), i=1
where the coefficients y0 and y1 are the solution of the linear system (17) with additional condition (18) being written down in block-matrix form (using regularization procedure similar to one described in 131):
A00 + D00 A01 + D01
A10 + D10 A11 + D11 O L0 L1 0
'y 0^ Y1 R )
b1 Vr/
Here, Apq are N x N square matrices; Dpq are diagonal matrices, p, q = 0, 1; I and O are vectors consist of ones and zeros, respectively; L0 and L1 are vectors consist of integrals from the corresponding basis functions along the curvilinear panels; yp = (y? , ..., yN)T are vectors of unknown coefficients; R is regularization variable; bp are right-hand side vectors:
A?/ = J tf(s) ds( J Q(s,a)^q(a) dA D?q = -1 J (s) ds,
Ki Vj ' Ki
L? = J ^P(s) ds = J f (s)^p(s) ds,
Ki Ki
i, j = 1, ..., N, p, q = 0, 1.
Basis functions are orthogonal, then
D01 = D 10
0
The length of the i-th panel ^ = Sj — Sj_i is considered to be small, i.e., we write down the approximate formulae for the matrix coefficients as series expansions with respect to ^ and take into account only the terms that are proportional to )k, where k ^ 3. The parameter value s* = Sj -1 corresponds to the center of the i-th panel.
For the diagonal coefficients of the matrices we obtain the following exact expressions:
D00 = — 2, D01 = D10 = 0, D11 = — A. (23)
The tangent unit vector t(s), the principal normal unit vector v(s) and the curvature
) |r'(s) x r''(s)|
K(S)= Ms)|>
of a planar curve are related by the Frenet — Serret formulae:
dT(S)= K*(s) v(s), ^VM = —K*(s) t(s). ds ds
But in the framework of this paper, the usage of outer normal unit vector n(s) and value
V(s) x r"(s)) ■ k
|r'(s)|3
which we call "signed curvature", is more convenient. Using the formulae similar to the Frenet — Serret formulae:
dT (s) f ^ f v dn(s) f \ f \
— =—K(s) n(s), — =K(s) T (s),
we can write down the expression for the derivatives of the radius-vector with respect to the natural parameter (arc length):
dr(s) , , d2r(s) . . d3r(s) .. . . . 2. . . . ^^
= t(s), = -k(s) n(s), = -k'(s) n(s) - K2(s) t(s). (24)
Numerical computation of the coefficients that form Aj matrices is not a trivial problem, because for the diagonal coefficients (i = j) and the coefficients which refer to the neighboring panels (|i — j | = 1), the denominator of the integrand is equal to zero at |r — £| ^ 0.
However, it is possible to obtain approximate analytical formulae for these coefficients: for the diagonal terms Aff they are the following:
x2 ¿3 r3
= Ki + O(i4), A0»1 = ï|4-Ki + O(i4), Ai = A-Ki + O(*4), A1/ = O(^4). (25)
Here, Ki is the curvature of the curve at the center of the i-th panel; the prime means the derivative with respect to the arc length.
Now we consider the coefficients which calculation requires integration over the different panels (i = j). We use the following additional designations hereinafter: dij = ri — Vj is the vector which connects the centers of the corresponding panels; ti, t j, ni, nj are tangent and outer normal unit vectors at the centers of the i-th and j-th panels, respectively.
The components of the matrix Aj (for i = j) are calculated according to the following rules: if the distance between the centers of the i-th and j-th panel is rather small, e.g.
|dij| ^ 2max{j,
the approximate formula should be used:
Aj - <dij ■ ni) + (((6 — — Kj (dij ■ nj ))(dij ■ ni) +
+ (Kj (dij ■ ti) — 4(tj ■ ni^ (dij ■ tj) j j + + ((6 — 8(dd i) — Ki (dij ■ (Kidij — 3ni))^ (dij ■ ni) —
— (dij ■ (3 KiT i — Ki dij )) (dij ■ t i^ b2i j . (26) For the panels placed far one from the other, the simplified formula is suitable:
Aj - (dij ^ ni>. (27)
Note that instead of (27) the value can be used which was derived for the rectilinear panels (as it was mentioned above, in this case the integral is computed analytically [25]).
For the components A0,1 and A1,0 we obtain
Aj ~ 24j 1V(Ti nj) +2 j 1, (28)
S2S, (d T W K 2(dij ■ ni) \ (29)
All the coefficients All, including diagonal, can be neglected.
Recall that the formulae (26) and (28)-(29) are correct not only for smooth but also for non-smooth airfoils, i.e., they can be applied for the neighboring panels, when there is the angle point of the airfoil at the common point of these panels.
However, in case of smooth airfoil, in practice the following approximate formulae, which are much simpler in comparison with (26) and (28)-(29), can be used for the matrix coefficients which calculation requires integration over the neighboring panels (|i — j | = 1).
/100 Mj/ Sj — 2Si „ j n /r4\ /101 _SiSj_ „J,n fX4\
A0 = ( Kij ±-6-/ + O(Sj), A°j = 144^Kij + O(Sj),
Alj0 = 72, Kij + Oj All = Oj
Here, Kij denotes the curvature of the airfoil at the common point of the neighboring panels; the prime, as earlier, means the derivative with respect to the arc length; sign "+" in the formula for A00 is used for j = i + 1, and "—" for j = i — 1.
We also note that in case of the airfoil with sharp edges, the formulae for rectilinear panels derived in [25] can be suitable for the matrix coefficients corresponding to the panels which are adjacent to the angle points. Those formulae doesn't permit to take into account the curvilinearity of the panels, however, they are exact for rectilinear panels.
In order to calculate the right-hand side (for the considered particular case when there is no vorticity in the flow domain f (s) = — V— ■ t(s)), the following approximate formulae are used:
b0 = —(V- ■ Ti)Si + ^(V- ■ (Kn + K2Ti))Si3 + O(Si4), (30)
bl = V- ■ niKiSi2 + O(Si4). (31)
Here, all the values are calculated at the centers of the corresponding panels; t and n are, as earlier, the tangent and outer normal unit vectors, respectively.
The results are shown in Figure 10 and 11, again only for the regions in the neighborhood of the left and right quadrant points of the ellipse.
It is seen that the obtained results are in good agreement with the most accurate solution obtained by numerical calculation of all the integrals along the exact curvilinear shape of the panels.
(a) (6)
Figure 10: Piecewise-constant(a) and piecewise-linear(ô) numerical solution in the neighborhood of the left quadrant point of the ellipse for the curvilinear panels at non-uniform ellipse discretization and semi-analytical coefficients computation
y(s) y(s)
(a) (6)
Figure 11: Piecewise-constant (a) and piecewise-linear (b) numerical solution in the neighborhood of the right quadrant point of the ellipse for the curvilinear panels at non-uniform ellipse discretization and semi-analytical coefficients computation
4.3. Approximate right-hand side integrals computation
In order to present the approximation technique completely, we derive the approximate expressions for all the terms in the right-hand side of the equation (17) (as earlier, only in the framework of the T-scheme, with above introduced basis and projection functions ^p = ^p, i = 1, ..., N, p = 0, 1), considering the vorticity in the flow domain to be given as a set of vortex elements,
N
)= E TwS(r — Pw).
w=1
In the most general case, the right-hand side, according to (9) and (17), instead of (30)—(31) has the following form:
b? = J f (s)tf(s) ds =
Ki
N N r
d ■ T(s)^(s) ds - E rj (Gw (s) ■ n(s)) <^(s) ds -
y w=1 /
VK2r,t) - V^ ■ T(s)^f(s) ds - > : r,„ / ( Gw(s) ■ n(
Ki w=! Ki
N m
-EE I tf (s) ds ( / f(G(s, a) ■ n(s)) (Yatt j + (G(s, a) ■ t(s)) (qatt j ^(a) da'
j=! q=0Ki y, v J '
or, the same,
N N m N m
bp = (bv)P - E rw(zw)P - EE Aj(yattj - EESj(qatt j, (32)
w=! j=!q=0 j=!q=0
where
(by)P = / - V«) ■ t(s)^P(s) ds,
Ki
(Zw)p = J (Gw(s) ■ n(s^ ^p(s) ds,
Ki
Sj = / tf (s) ds ( J(G(s, a) ■ t(s))¥j(a) da) ,
Ki K,
expressions for A j are the same as earlier.
For the coefficients (by)P the following approximate formulae are obtained:
(bv)0 = 2((VK,i - 2V^) ■ t>) + 4_((VK,i - 2V^) ■ Ti)" +
2 • - ■ ■ ■ 4g
£2 , r2 ✓ X
(by)! = 24((VK,i - 2V^) ■ Ti) + O(i4) = 24[(VK)i ■ Ti) - **((VK,i - 2V«) ■ n*)j + O(#).
Here, V K,j is the airfoil surface line velocity at the center of the i-th panel; t j and v j are the tangent and principal normal unit vectors at the center of the i-th panel, respectively; prime denotes the derivative with respect to the arc length.
Now we write down the approximate formulae for the coefficients . The technique of their acquisition is the same in principle as for the coefficients, however, the transformations is now more complicated since the integrand becomes singular even for smooth airfoil surface line segments, so the integral should be calculated as its principal value in Cauchy sense. For the diagonal (i = j) terms S?, the approximate formulae are the following:
1 X2
S00 = Sij1 = S0j1 = -S1° = - + 288n ^ +
When considering the coefficients which calculation requires integration over the different panels (i = j), we, as earlier, use the designations dj = r - r j for the vector which connects the centers
of the i-th and j -th panels. The coefficients Sj (i = j) are calculated according to the following rules: if the distance between the centers of the i-th and j -th panel is rather small, e.g.
| ^ 2max{j,
the approximate formulae should be used:
sj(dij ■ T o+-j (dij ■ T M j)2 -1 k+( jf -1 kl +
ij 2n|dij|2V ij ^ 2n|dij |4V ij ^V 3|dij|2 4) i ' V 3|dij|2 4) j . , \ . Wjr^Jr,! n x2 (dij ■ nj)kjj '
+ K (Ti ^ nJ) + 8^ (dij ^ Ti) ^ --3-,
S S3 S
i j |4 (tj ■ ni)(dij ■ n) - 48niidj p((dij ■ ni)Ki + (dij ■T¿)k2) ; (33)
12n|dij|4' j ^ ij j 48n|dij|
S01 _ ^ f (dij ■ Ti)(dij ■ Tj) (ti ■ T j A +
ij ~ 12n|dij |2^ |dij|2 - ' +
+ 8& (a|j((dij ■ Ti)3(dij ■ Tj) + (dij ■ ni)3(dij ■ nj^ - +
+ Jjû Ki f + (dij ■ nj ) ((dij ■ T ■T j ) - ^ )) +
+ 57*^, - 2j j ■ (Ki ni + K2Ti); (34)
S10 ~ ^ A _ (dij ■ ti)2 (dij ■ ni) K \ _
ij
12n|dij |2\ 2 |dij|2 2
"ij | \ | "f?
$ j f 7 (ti ■ Tj)2 (dij ■ ti)2 (dij ■ Tj)2
— - ■ - — - —
+ v ; , ,, +
4n|dij |4V 24 4 3|dij |2 3|dij |2 ,
I (d.. . T ,)(d.. . n N f (dij ■ T j )(dij ■ ni) (
+ 12n|dij|6 (dij Tj)(dij ni^ |dij |2 - (Ti nj V -
Si2 j f (ti ■ nj ) 1 A
-(KiTj + KjTi) — (Kini - Kjnj M ■ dij -
48n|dij|4\ 3 ^ ' j 1 j " 2
¿i^Kkj / (ti ■ tj) (dij ■ ni)(dij ■ nj)'
— - ■ - —
288n|dij |2\ 2 |dij |2
(35)
S11 _ S2S2 (d T (dij ■ ti)2 (dij ■ ni) A +
- 36,(dij ■Tj) u — r~ Kv +
S,2S2 . / (d, • ni) xA „ „
+ 72j <T i'nj H^jF — 4> (36)
For the panels placed far one from the other, the simplified formulae are suitable:
c00 _ Q01 ^ ^ Adij' ■ T i)(dij ■ T jN (t i ■ T j y
Sij ^„r12 (dij ^ TiN' Sij
2n|dijij ^ ij 12n|dij |dij|2
SA (dij ■ Ti)2 (dij ■ ni) \ S0 Sj - ^^ i2 - ^^ — KV ' Sj - 0-
As for coefficients, the formulae (33)-(36) remain correct for non-smooth airfoils, i.e., they can be applied for the panels when there is the angle point of the airfoil between them.
For smooth parts of the airfoil the following approximate formulae, which are much simpler, but more accurate in comparison to (33)-(36), should be used for the matrix coefficients calculation that requires integration over the neighboring panels (|i — j| = 1):
S
00 ij
J— * in \ 2n
* +
— *L in 2n
) 2n V
* + j , № + j + °(*4),
+
48n
Soi = _ + + ) lnf * + A +
281 ^ +
S10 = j _ (* + ) i^ * +
ij
S
ii
ij
±
4n*i * + *j
24n
inf * +
*2*
288n "ij
i 4 + °(*4)
24n*
lnf * +
24n*
+ °(*4).
Here, Kij denotes the curvature of the airfoil at the common point of the neighboring panels; sign "+" in the formulae for Sij0 and Si1 are used for j = i + 1, and "—" for j = i — 1.
We note again that in case of the airfoil with sharp edges, the formulae for rectilinear panels derived in [25] can be suitable for the matrix coefficients corresponding to the panels which are adjacent to the angle points. Those formulae don't permit to take into account the curvilinearity of the panels, however, they are exact for rectilinear panels.
In order to take into account the influence of the vortex elements in the flow, i.e., to calculate coefficients (zw)p in the formula (32), two cases should be considered. Let us denote the vector which connects the position of the w-th vortex element with the center of the i-th panel as hiw =
In case of small
5i
I hi
ratio, i.e., when the vortex element is rather far from the panel, the
following approximate formulae are suitable:
0 _ (hiw ' ni) *
w)i = 2n|hiw| | ' 2n\ V4
(Zw )
+ (hiw ■ ni)2 K2|hiw|2\ (hiw ■ ni) +
31 hiw 12
24
Ki 1 hiw 1 Ki ( hiw ■ tih (hiw ■ ti) + Ki|hiw^f *
24
4
(zw )1 =
|hiw | 3
|hi
8
. |hi
+ °
(hiw ■ T i) ( Ki | hiw | (hiw ■ ni) A ( *i
12n|hi
2
I hi
| hiw |
| hiw |
2 / *
+
hiw
If the vortex element is placed close to the panel, it is hardly possible to derive rather simple and at the same time rather accurate approximate formulae, taking into account the curvilinearity of the panel. However, in this case the panel can be replaced with the arc of the osculating circle of radius Ri = — having the same length as the original panel.
Under this assumption, it is possible to use the exact solution [26] and write down the following formulae:
(zw )0
( Xi +
( Xi
2
4
4
where
a,
,TI0 ß - ft* , 1 f f f f ß - ft
^(ß+ 2n — 2 + a,
- n sign a,wn (ß - ft*)
is the continuous branch of the antiderivative for the corresponding integrand with respect to the arc length; x- is the argument of the (r — c) vector (the angle between this vector and horizontal
line, Figure 12); vector c- = r--- denotes the center of the osculating circle; - and k-, as
earlier, are the length of the i-th panel and the signed curvature at its center; k* = | k | ^ 0 is the curvature; is the argument of the (pw — c-) vector, which value can be corrected by ±2n in order to provide the condition satisfaction
—K* A. K*
—n<P — < n for Xi — At < P < x- + -y-;
aiw is dimensionless parameter denoting how close the vortex element is placed to the arc length,
- 1= |pw - C1 -
R,
n (x) is the Heaviside function, which is equal to 1 at x > 0 and equal to 0 at x < 0
Figure 12: Vortex element position, airfoil curvilinear panel and its approximation with osculating circle
The developed approach doesn't permit to calculate (zw)1 coefficient. To be more precise, it is possible to integrate the corresponding function along the arc of the osculating circle, but the result is expressed through the polylogarithmic special functions, which are not suitable for practical purposes. However, considering the above introduced parameter aiw to be small, it is possible to expand the polylogarithms (as well as other rather complicated terms) into power series with respect to aiw and write down the following approximate formula for (zw)1:
(zw)i ~ ^¿w ( Xi + 2 ) ( X 2 /
where (ß) is the continuous branch of the antiderivative for the corresponding integrand with respect to the arc length of the osculating circle:
(ß) =
1
k*
(ß - ft*):
+
+ I (P — ) — 2nn(aiW) — 2 arctan
P a2w (P ^¿w)
Xi) + 5aiw — 18
+ (P — tfiw)
2
+ aiw (2 aiw )
2
18 aw + (p—^¿w)2
Conclusion
The problem of the accuracy improving is considered for vortex method. The general Galerkin-type approach is considered for the numerical solution of the boundary integral equation. It is shown that the airfoil surface line representation as a polygon consisting of rectilinear panels can lead to incorrect behavior of the numerical solution of boundary integral equation with respect to unknown vortex sheet intensity, especially in case of essentially different lengths of the neighboring panels. However, in this case the exact analytical formulae are known for the matrix coefficients which approximate the integral equation. The approach is developed which makes it possible to take into account explicitly the curvature of the surface line. The approximate analytical formulae are derived for the coefficients of the linear system which are expressed through the integrals over the curvilinear panels. Numerical example is considered for the elliptic airfoil in the incident flow at its essentially non-uniform discretization. The obtained result are very close to the most accurate one, which can be obtained by numerical calculation (too complicated) of all the integrals along the exact curvilinear shape of the panels.
This study was partially funded by the Russian Foundation for Basic Research according to the research project No. 18-31-20051.
1. Cottet G.-H., Koumoutsakos P.D. Vortex methods: Theory and practice. Camb.; N.Y.: Camb. Univ. Press, 2000. 313 p.
2. Lewis R.I. Vortex element methods for fluid dynamic analysis of engineering systems. Camb.; N.Y.: Camb. Univ. Press, 2005. 566 p.
3. Lifanov I.K. Singular integral equations and discrete vortices. Utrecht: VSP, 1996. 475 p.
4. Lifanov I.K., Poltavskii L.N., Vainikko G.M. Hypersingular integral equations and their applications. Boca Raton; L.: Chapman & Hall: CRC Press, 2004. 396 p.
5. Ostrikov N.N., Zhmulin E.M. Vortex dynamics of viscous fluid flows: Pt. 1: Two-dimentional flows. J. of Fluid Mechanics, 1994, 276, pp. 81-111. DOI: D0IURL10.1017/S0022112094002478
6. Dynnikova G.Ya. The Lagrangian approach to solving the time-dependent Navier - Stokes equations. Doklady Physics, 2004, vol. 49, no. 11, pp. 648-652. DOI: 10.1134/1.1831530
Acknowledgments
References
7. Guvernyuk S.V., Dynnikova G.Ya. Modeling the flow past an oscillating airfoil by the method of viscous vortex domains. Fluid Dynamics, 2007, vol.42, no. 1, pp. 1-11. DOI: 10.1134/S0015462807010012
8. Dynnikov Ya.A., Dynnikova G.Ya. Application of viscous vortex domains method for solving flow-structure problems. ECCOMAS Thematic conf. on multibody dynamics (Zagreb, Croatia, July 1-4, 2013): Proc. Zagreb, 2013. Pp. 877-882. DOI: 10.13140/2.1.2113.1207
9. Wu J.C., Thompson J.F. Numerical solutions of time-dependent incompressible Navier-Stokes equations using an integro-differential formulation. Computers & Fluids, 1973, vol. 1, no. 2, pp. 197-215. DOI: 10.1016/0045-7930(73)90018-2
10. Bykhovskii E.B., Smirnov N.V. Orthogonal decomposition of the space of vector functions square-summable on a given domain, and the operators of vector analysis. Trudy Matem-aticheskogo Insituta im. Steklova RAS [Proc. of Steklov Math. Inst.], 1960, vol. 59, pp. 5-36. Available at: http://mi.mathnet.ru/eng/tm/v59/p5, accessed 3.12.2018 (in Russian).
11. Kempka S.N., Glass M.W., Peery J.S., Strickland J.H., Ingber M.S. Accuracy considerations for implementing velocity boundary conditions in vorticity formulations. SANDIA Report, 1996, no. SAND96-0583 UC-700. DOI: 10.2172/242701
12. Rosenhead L. The formation of vortices from a surface of discontinuity. Proc. of the Royal Soc. of London. Ser. A: Math., Physical and Engineering Sciences, 1931, vol. 134, no. 823, pp. 170-192. DOI: 10.1098/rspa.1931.0189
13. Lighthill M.J. Introduction. Boundary layer theory. Laminar boundary layers / Ed. by L. Rosenhead. Oxf.: Clarendon Press, 1963. Ch. 2. Pp. 46-113.
14. Muskhelishvili N.I. Singular integral equations. Groningen: P. Noordhoff, 1953. 447 p. (Russ. ed.: Muskhelishvili N.I. Singuliarnye integral'nye uravneniia. Moscow: Gostekhizdat Publ., 1946. 448 p.).
15. Dynnikova G.Ya., Andronov P.R. Expressions of force and moment exerted on a body in a viscous flow via the flux of vorticity generated on its surface. Eur. J. of Mechanics { B/Fluids, 2018, vol. 72, pp. 293-300. DOI: 10.1016/j.euromechflu.2018.06.002
16. Dynnikova G Ya. The integral formula for pressure field in the nonstationary barotropic flows of viscous fluid. J. of Mathematical Fluid Mechanics, 2014, vol. 16, no. 1, pp. 145-162. DOI: 10.1007/s00021-013-0148-z
17. Yoshifumi Ogami, Teruaki Akamatsu. Viscous flow simulation using the discrete vortex model — the diffusion velocity method. Computers & Fluids, 1991, vol. 19, no. 3-4, pp. 433441. DOI: 10.1016/0045-7930(91)90068-S
18. Dynnikova G.Ya. Vortex motion in two-dimensional viscous fluid flows. Fluid Dynamics, 2003, vol. 38, no. 5, pp. 670-678. DOI: 10.1023/B:FLUI.0000007829.78673.01
19. Moreva V.S., Marchevsky I.K. Vortex element method for 2D flow simulation with tangent velocity components on airfoil surface. 6th Eur. congress on computational methods in appliedsci-
ences and engineering: ECCOMAS 2012 (Vienna, Austria, September 10-14, 2012): Full papers. Vienna, 2012. Pp. 5952-5965. Available at: http://eccomas.org/cvdata/cntr1/spc7/dtos/, img/mdia/ECCOMAS-2012-e-book-Title-Content.pdf, accessed 3.12.2018.
20. Kuzmina K.S., Marchevsky I.K., Milani D., Ryatina E.P. Accuracy comparison of different approaches for vortex sheet discretization on the airfoil in vortex particles method. Particles 2017: 5th Intern. conf. on particle-based methods — Fundamentals and applications (Hannover, Germany, September 26-28, 2017): Proc. 2017. Pp. 691-702. Available at: http://www.eccomas.org/cvdata/cntr1/spc10/dtos/img/mdia/particles2017-ebook.pdf, accessed 3.12.2018.
21. Kuzmina K.S., Marchevskii I.K., Moreva V.S. Vortex sheet intensity computation in incompressible flow simulation around an airfoil by using vortex methods. Mathematical Models and Computer Simulations, 2018, vol. 10, no. 3, pp. 276-287. DOI: 10.1134/S2070048218030092
22. Kuz'mina K.S., Marchevskii I.K., Moreva V.S., Ryatina E.P. Numerical scheme of the second order of accuracy for vortex methods for incompressible flow simulation around airfoils. Russian Aeronautics, 2017, vol. 60, no. 3, pp. 398-405. DOI: 10.3103/S1068799816030114
23. Katz J., Plotkin A. Low-speed aerodynamics: From wing theory to panel methods. N.Y.: McGraw-Hill, 1991. 632 p.
24. Cockburn B., Shu C.-W. Runge — Kutta discontinuous Galerkin methods for convection-dominated problems. J. of Scientific Computing, 2001, vol.16, no. 3, pp. 173-261. DOI: 10.1023/A:1012873910884
25. KuzminaK.S., Marchevsky I.K., RyatinaE.P. Exact analytical formulae for linearly distributed vortex and source sheets influence computation in 2D vortex methods. Journal of Physics: Conference Series, 2017, vol. 918, no. 1, art. no. 012013. DOI: 10.1088/1742-6596/918/1/012013
26. Kuzmina K.S., Marchevsky I.K., Ryatina E.P. Exact solutions of boundary integral equation arising in vortex methods for incompressible flow simulation around elliptical and Zhukovsky airfoils. Journal of Physics: Conference Series, 2019 (in press).
Ссылка на статью:
// Математика и математическое моделирование. 2018. №6. С. 22-51.
Б01: 10.24108/шаШш.0618.0000158
Представлена в редакцию: 19.11.2018 © НП <<НЭИКОН>>
Математика U Математическое
моделирование
Сетевое научное издание http://mathmelpub.ru ISSN 2412-5911
УДК 004.3+519.6
Уточненный алгоритм аппроксимации граничного интегрального уравнения в вихревых методах моделирования обтекания профилей с криволинейной границей
Кузьмина К. С.1'2, Марчевский И. К.1'2, Солдатова И. А.1
1МГТУ им. Н.Э. Баумана, Москва, Россия 2Институт системного программирования имени В.П. Иванникова РАН, Москва
Ключевые слова: вихревой метод; граничное интегральное уравнение; плоское течение; криволинейная панель; численное интегрирование; неравномерная дискретизация
■к
Рассматривается проблема повышения точности лагранжевых вихревых методов моделирования обтекания профилей плоскопараллельным потоком вязкой несжимаемой среды. Определяющим фактором является точность моделирования процесса генерации завихренности на границе профиля. В рассмотренной математической модели слой завихренности, генерируемой за шаг расчета по времени, представляется тонким вихревым слоем, расположенным на границе области течения - на поверхности обтекаемого профиля. Его интенсивность может быть найдена из граничного условия прилипания, выраженного в форме векторного граничного интегрального уравнения, при этом достаточно обеспечить его удовлетворение лишь для одной из компонент: нормальной либо касательной к границе профиля. Используемый математический аппарат основан на свойствах обобщенного разложения Гельмгольца для соленоидальных векторных полей. Получающиеся в результате скалярные интегральные уравнения имеют существенно различные свойства: в первом случае уравнение является сингулярным и интеграл в нем понимается в смысле главного значения по Коши, во втором - ядро уравнения ограничено для гладких участков границы. В работе рассмотрены оба подхода.
Для численного решения граничных интегральных уравнений в вихревых методах обычно используют так называемые панельные методы, в соответствии с которыми исходная кривая, задающая границу профиля, заменяется панелями, для которых записывают дискретный аналог исходного интегрального уравнения. В работе предложен подход к аппроксимации граничного интегрального уравнения, основанный на идеях разрывного метода Га-
леркина. При этом в качестве базисных и проекционных функций могут быть использованы дельта-функции Дирака, кусочно-константные и кусочно-линейные функции. Показано, что дискретизация профиля прямолинейными панелями может приводить к существенной погрешности приближенного решения, в особенности в случае существенно различных длин соседних панелей. Однако для такой дискретизации получены точные аналитические формулы для вычисления интегралов, выражающих коэффициенты системы линейных уравнений, аппроксимирующей интегральное уравнение в соответствии с подходом Галеркина.
Предложен подход, позволяющий явно учесть кривизну границы профиля. Для коэффициентов системы линейных алгебраических уравнений, представляющей дискретный аналог интегрального уравнения, получены приближенные аналитические выражения в виде степенных разложений до слагаемых, пропорциональных третьей степени малого параметра — длины панели. Рассмотрен пример численного моделирования обтекания эллиптического профиля при существенно неравномерной дискретизации границы профиля. Полученный алгоритм имеет приемлемую вычислительную сложность и при этом позволяет получать верное качественно и более точное количественно решение по сравнению с ранее известными подходами.
Список литературы
1. Cottet G.-H., Koumoutsakos P.D. Vortex methods: Theory and practice. Camb.; N.Y.: Camb. Univ. Press, 2000. 313 p.
2. Lewis R.I. Vortex element methods for fluid dynamic analysis of engineering systems. Camb.; N.Y.: Camb. Univ. Press, 2005. 566 p.
3. Lifanov I.K. Singular integral equations and discrete vortices. Utrecht: VSP, 1996. 475 p.
4. Lifanov I.K., Poltavskii L.N., Vainikko G.M. Hypersingular integral equations and their applications. Boca Raton; L.: Chapman & Hall: CRC Press, 2004. 396 p.
5. Ostrikov N.N., Zhmulin E.M. Vortex dynamics of viscous fluid flows: Pt.1: Two-dimensional flows // J. of Fluid Mechanics. 1994. Vol. 276. Pp. 81-111. DOI: 10.1017/S0022112094002478
6. Дынникова Г.Я. Лагранжев подход к решению нестационарных уравнений Навье — Стокса // Докл. Российской акад. наук (РАН). 2004. Т. 399, № 1. С. 42-46.
7. Гувернюк С.В., Дынникова Г.Я. Моделирование обтекания колеблющегося профиля методом вязких вихревых доменов // Изв. РАН. Механика жидкости и газа. 2007. №1. С. 3-14.
8. Dynnikov Ya.A., Dynnikova G.Ya. Application of viscous vortex domains method for solving flow-structure problems // ECCOMAS Thematic conf. on multibody dynamics (Zagreb, Croatia, July 1-4, 2013): Proc. Zagreb, 2013. Pp. 877-882. DOI: 10.13140/2.1.2113.1207
9. Wu J.C., Thompson J.F. Numerical solutions of time-dependent incompressible Navier-Stokes equations using an integro-differential formulation // Computers & Fluids. 1973. Vol. 1, no. 2. Pp. 197-215. DOI: 10.1016/0045-7930(73)90018-2
10. Быховский Э.Б., Смирнов H.B. Об ортогональном разложении пространства вектор-функций, квадратично суммируемых по заданной области и операторах векторного анализа // Тр. МИАН СССР. 1960. Т. 59. С. 5-36. Режим доступа: http://mi.mathnet.ru/ eng/book1108 (дата обращения 3.12.2018).
11. Kempka S.N., Glass M.W., Peery J.S., Strickland J.H., Ingber M.S. Accuracy considerations for implementing velocity boundary conditions in vorticity formulations // SANDIA Report. 1996. No. SAND96-0583 UC-700. DOI: 10.2172/242701
12. Rosenhead L. The formation of vortices from a surface of discontinuity // Proc. of the Royal Soc. of London. Ser. A: Math., Physical and Engineering Sciences. 1931. Vol. 134, no. 823. Pp. 170-192. DOI: 10.1098/rspa.1931.0189
13. Lighthill M.J. Introduction. Boundary layer theory // Laminar boundary layers / Ed. by L. Rosenhead. Oxf.: Clarendon Press, 1963. Ch. 2. Pp. 46-113.
14. Мусхелишвили Н.И. Сингулярные интегральные уравнения. М.; Л.: Гостехиздат, 1946. 448 с. [MuskhelishviliN.I. Singular integral equations. Groningen: P.Noordhoff, 1953.447p.].
15. Dynnikova G.Ya., Andronov P.R. Expressions of force and moment exerted on a body in a viscous flow via the flux of vorticity generated on its surface // Eur. J. of Mechanics - B/Fluids. 2018. Vol. 72. Pp. 293-300. DOI: 10.1016/j.euromechflu.2018.06.002
16. Dynnikova G.Ya. The integral formula for pressure field in the nonstationary barotropic flows of viscous fluid//J. of Mathematical Fluid Mechanics. 2014. Vol. 16, no. 1. Pp. 145-162. DOI: 10.1007/s00021-013-0148-z
17. Yoshifumi Ogami, Teruaki Akamatsu. Viscous flow simulation using the discrete vortex model — the diffusion velocity method// Computers & Fluids. 1991. Vol. 19, no. 3-4. Pp. 433441. DOI: 10.1016/0045-7930(91)90068-S
18. Дынникова Г.Я. Движение вихрей в двумерных течениях вязкой жидкости // Изв. РАН. Механика жидкости и газа. 2003. № 5. С. 11-19.
19. Moreva V.S., Marchevsky I.K. Vortex element method for 2D flow simulation with tangent velocity components on airfoil surface // 6th Eur. congress on computational methods in applied sciences and engineering: ECCOMAS 2012 (Vienna, Austria, September 10-14, 2012): Full papers. Vienna, 2012. Pp. 5952-5965. Режим доступа: http://eccomas.cimne.com/ cvdata/cntr1/spc7/dtos/img/mdia/ECCOMAS-2Q12-e-book.pdf (дата обращения 3.12.2018).
20. Kuzmina K.S., Marchevsky I.K., Milani D., Ryatina E.P. Accuracy comparison of different approaches for vortex sheet discretization on the airfoil in vortex particles method // Particles 2017: 5th Intern. conf. on particle-based methods — Fundamentals and applications
(Hannover, Germany, September 26-28, 2017): Proc. 2017. Pp. 691-702. Режим доступа: http://www.eccomas.org/cvdata/cntr1/ spc10/dtos/img/mdia/particles2017-ebook.pdf (дата обращения 3.12.2018).
21. Кузьмина К.С., Марчевский И.К., Морева B.C. Определение интенсивности вихревого слоя при моделировании вихревыми методами обтекания профиля потоком несжимаемой среды // Математическое моделирование. 2017. Т. 29, № 10. С. 20-34.
22. Кузьмина К.С., Марчевский И.К., Морева B.C., Рятина Е.П. Расчетная схема вихревых методов второго порядка точности для моделирования обтекания профилей несжимаемым потоком // Изв. высших учебных заведений. Авиационная техника. 2017. №3. С. 73-80.
23. Katz J., Plotkin A. Low-speed aerodynamics: From wing theory to panel methods. N.Y.: McGraw-Hill, 1991. 632 p.
24. Cockburn B., Shu C.-W. Runge — Kutta discontinuous Galerkin methods for convection-dominated problems // J. of Scientific Computing. 2001. Vol. 16, no. 3. Pp. 173-261. DOI: 10.1023/A:1012873910884
25. Kuzmina K.S., Marchevsky I.K., Ryatina E.P. Exact analytical formulae for linearly distributed vortex and source sheets influence computation in 2D vortex methods // Journal of Physics: Conference Series. 2017. Vol.918, no. 1. Art. no. 012013. DOI: 10.1088/1742-6596/918/1/012013
26. Kuzmina K.S., Marchevsky I.K., Ryatina E.P. Exact solutions of boundary integral equation arising in vortex methods for incompressible flow simulation around elliptical and Zhukovsky airfoils // Journal of Physics: Conference Series. 2019 (в печати).
