On one numerical method of integrating the dynamical equations of projectile planar flight affected by wind Текст научной статьи по специальности «Физика»

UDC 531.55

On One Numerical Method of Integrating the Dynamical Equations of Projectile Planar Flight Affected by Wind

V. V. Chistyakov

Faculty of Engineering Yaroslavl' State Academy of Agriculture 58, Tutaevskoe highway, Yaroslavl', Russian Federation, 150042

Common way to integrate the dynamical equations of projectile planar motion introduces two Cartesian coordinates x(t) and y(t) and attack angle all depending on time t,

and three coupled ordinary differential equations (ODE) each nominally of Il-nd order. It leads to inevitable computational complexities and accuracy risks. The method proposed excludes the time variable and diminishes the number of functions to n = 2: the attack angle •d(b) and intercept a(b) of the tangent to the trajectory at the point with slope b = tan 9 with the 9 being the inclination angle. This approach based on Legendre transformation makes the integration more convenient and reliable in the studied case of quadratic in speed aerodynamic forces i.e. drag, lifting force, conservative and damping momenta and the wind affecting the flight. The effective dimensionality of new ODE system is diminished by 2 units and its transcendence is eliminated by simple substitution ij = sin •&. Also the method enables to obtain easily and reliably the projectile trajectories in conditions of tail-, head-or side wind. Investigated are main ranges of aerodynamic parameters at which takes place different behavior of the attack angle •& vs slope b including quasi-stabilization and aperiodic auto-oscillations. In addition, it was revealed non-monotonous behavior of speed with two minima while projectile descending if launched at the angles 90 close to 90°. The numerical method may implement into quality improvement of real combat or sporting projectiles such as arch arrow, lance, finned rocket etc.

Key words and phrases: projectile, lifting force, quadratic drag, conservative/damping momenta, attack angle, projective-dual variables, wind.

1. Introduction

In stationary conditions of no wind and when neglecting the Coriolis force the motion of finned symmetrical rigid projectile launched in the plane of its symmetry will stay planar up to its landing assumed as a rule at the start height. The tail-or head wind with constant velocity changes the projectile trajectory but conserves flight plane. As for side wind, it turns the plane around vertical direction on some angle which value depends on speeds' ratio, flight time and so on.

The dynamic system describing the flight includes n = 3 ODEs with two ones in Cartesian coordinates being of order k = 1 responding for mass-center (c. m.) motion and the third of order k = 2 describing the rotation of launched projectile in the flight plane around the c. m.

As a typical example of such dynamic system may be taken arch arrow [1] or lance, or finned combat projectile shot from the gun, or rocket moving freely after the short correcting jet impulse [2]. Also for the bullet shot from smoothbore gun. In all these and other analogous cases, the flight trajectories of c. m.'s and their characteristics are of large more interest than detailed flight process in real time. In addition, the projectile orientation in landing moment is not less importance. This is especially for combat projectiles in order to avoid ricochets and for effective target engagement.

Therefore the problem of integrating of ODEs of planar resistive motion especially in wind condition is actual both in mechanics of archery/lance and in exterior ballistics of finned projectiles. Moreover, the aerodynamically symmetrical projectile is a relatively simple experimental model for verifying qualitative and even quantitative conclusions of the theory of so-called variable dissipation systems [3] with its complex and various mathematical apparatus.

Received 6th June, 2014.

In this paper, an alternative way is developed how to determine the trajectory of projectile motion in vertical plane of its own symmetry also the attack angle behavior in assumption of quadratic law for all aerodynamic efforts regarded drag and lifting forces, conservative and damping momenta. Earlier the method was tested successfully for a heavy point [4,5] and gyro-stabilized axially symmetrical rigid projectile [6] and it is based on Legendre transition from Cartesian coordinates (x(t), y(t)) to dual projective variables. These ones are slope of the c. m. trajectory b = tg 0 and intercept a = y — bx or more precisely its variation.

The temporal characteristics considered are in second turn in this work as not being of direct interest but influencing some physical values, e.g. a spectrum of infrasonic waves emitted on different parts of the trajectory.

An additional advantage of dual-projective variables is their monotonic behavior and invariance relatively parallel shift of the coordinate axes.

As for the slope b = tg 0 this value is more preferable as natural angular measure in ballistics than the angle d itself because of more convenience of targeting by ratio of vertical and horizontal catheta than that of arc and radius.

2. Primary Equations with no Wind Affecting

It's necessary first to integrate the flight dynamical equations at no macroscopic air motion or the same in wind-connected frames assumed be inertial due to constant wind velocity U = (u, 0, w). The solution obtained should then be adapted for fixed frames connected with launching point and shot direction.

It is assumed that the projectile has some longitudinal axis I of symmetry along which is directed the initial velocity V0 = (V0 cos d0,V0 sin d0, 0) with d0 = 0(0) being the angle of throwing, and drag force is minimal for this direction too. In addition, the lift force is equal to zero for the direction though this is not obligatory and the body may be launched at attack angle of $0 = 0. The last takes place when taking into account tail- or headwind.

The projectile rotation considered is in non-inertial reference frames with its beginning in mass center point C so the gravitational and entrainment inertial momenta are equal to zero. Therefore, the dynamics of attack angle defined is only by aerodynamic torque Mc decomposing on conservative static part Ms and the dissipative damping momentum M^.

For the momenta above used are model formulae [7] used in exterior ballistics as follows Ms = —aV2 sin Ma = —5V2w where V is velocity of c. m., 6 for the attack angle, u for angular velocity, a and 5 both are the corresponding coefficients assumed in the model not depending on § and w.

As for the motion of mass center itself observed from fixed inertial reference frames Oxy (Fig. 1) it is defined along with gravity mg also by the aerodynamic force vector

Rad decomposing on the drag R = — mga0(1 + £ sin2 §)V2 and normal force N = mg^0V2 sin $ both quadratic in speed V. The finned tail makes it impossible for the

projectile to rotate about longitudinal axis I, so the Magnus force and momentum don't appear.

The parameter s describing relative increase of the drag R under changing of attack angle from 0 to 90° may be large enough, i.e. for arrow of sporting arch this may achieve some decades because of the fact that along with large increase of the frontal square S hugely worsens its streamlining.

In fixed reference frames Oxy the dynamic equations for mass center are as follows

/-, ■ 2 n\ 9%2 ■ n9%2 sin 0

x = — «o(1 + £ sin v)-- — 70 sin §-r—-,

cos 0 cos2 0

/-, ■ 2 n\ 9%2 sin 0 . . gx2 (1)

y = —«o(1 + £ sin -&)-+ 70 sin §-- — g,

cos2 0 cos 0

^ y = tan 0 ■ x. The system may be reduced to

/-, ■ 2 n\ 9%2 ■ n9%2 sin 0

x = —ao (1 + £ sin v)-- — 70 sin §-——,

cos 0 cos2 0

sin -dgx g cos2 6

6 = 70

cos d

In frames moving in translational way with c. m. rotation of the projectile is described by the ODE of Il-nd with respect to the orientation angle 0 = 6 + § between the axis I and horizon Ox

j d2(d + #) i2 i2 d(0 + ,0, = —^co^ sin * — °cos2 0 di , (2)

where Jcz is the central transversal momentum of inertia, i.e. with respect to horizontal axis Cz.

Thus, we receive the ODE system with three unknowns that is complicated enough. By solving it, we find time dependences for attack § and inclination angle d and for horizontal velocity x too. The nominal order of that system is equal k = 5 but in fact is 4.

3. New Equation for Mass Center

y

In Eqs (1)-(2) when taking into account the relations b = tan d = — and 0 = d + $

there are three independent spatial variables and argument t — the time past from the start.

Multiplying the first equation in (1) by b and then subtracting it from the second one it receives

y — bx = —g + 70 sin ■dgx2(1 + b2)3/2. Given the fact that y = bx it means

bx = —g + 7 (b)gx2(1 + b2)3/2.

As far as x = — da(^ [6] and hence x = — <d2rb > 0 then do do2

& (—d2i *)=—^+^o sin ^ (—d2i &)2 (i+^

(3)

and

i

b + ™ sin 0(1 + b2)3/2 • (4)

The negative sign is ascribed there when normal gravity projection exceeds the lifting force mg cos 6 > N, otherwise curvature of trajectory is negative and projectile flies convex trajectory.

Mass center velocity determined is as

v (b) = £¿(1 + &2)1/2 = ^ 9(1 +b2)<b

db Y 1 + 70 sin ïïga'lb(1 + b2)3/2

(

9(1 + b2) (5)

^ + 7o sin + 62)3/2 ' ( )

For the square of loss of the ^-inclination rate it receives

¿2 =_g_= ___1_ (6)

_l_ oir, Sin" \2H_Lh2\3/2 („» \2 „o/i , L2V3/2' W

bb

a'

+ 970 sin 0(a'b'b)2(1 + b2)3/2 (a'bb)2 ^ + g^ sin 0(1 + 62)3/2'

For the derivative of II-nd order respectively

9 Uib + 297o sin ^a^a + b2)3/2 + ^(ag^)2d(7° sin^+fe2)3/2))

b =----2--. (7)

2«; + 97o sin 0(<6)2(1 + 62)3/2)2 Horizontal acceleration of c. m. is determined in new variables as

( d2aA d3a- 2 d2av

= — d*"2 — <8> Then the first of Eqs (1) is written after substitution (4), (6), (7) as

_ n"> _I__

bbba'lb + 70 sin ^(al)2(1 + b2)3/2

/

abb

9 (a'»b + 270 sin Vg^bb (1 + b2)3/2 + ffW^ sin^+b2)3/2) ) ^ 2(ab + 7o sin «b)2(1 + 62)3/2)2 ^

= — (ao(1 + £ sin2 0) + 706 sintf)g(—a'{>b)2x

a'bb + 7o sin Vg«b)2(1 + b2)3/2 After simple rearrangement, it takes the form

+ b2. (9)

a*hh = 2g{a(b) + by (b))abb (4" + 9l(b)(1 + b2 )3/4 Vl+b2+

\ab b J

(al)2 w ,w/bbV«:

d( 7 (6)(1 + b2)3/2) + 3--. (10)

The a(b) and 7(b) are projectile drag and lift force coefficients respectively in very general case but here these a(b) = ao (1+£ sin2 §(b)) and 7(6) = 70 sin §(b). Therefore, the equation is

obbb , ^^^„^„/^ ± + o7(6)(1 + 62)3/M V1 +

= 2g(ao(1 + e2m) + lob sin ti(b))<b{ 4" + 91 (&)(! + 62)3/^ + b2+

\ab b /

Kb)2 Va

d(7o sin tf(&)(1 + b2)3/2) / + 9—-^-^. (11)

Its solution would describe parametrically the trajectory and flight time plot through following inverse transition formulas [4,5]:

b b f d2a(b')A [ d2a(b')A

x(b) — xo = - J db , V(b) — yo = - J db ,

b f (12)

t(b) = — ^ ,--=, 6o = tan do.

bo V a'b',b, + 7o sin <d(V)g(a'l,b,)2(1 + b'2)3/2

The parameter b above gets useful geometrical meaning, namely inclination of vector V to the horizon. However, independent integration of the equation is impossible since it contains unknown function ^(b) = sin $ = sin (0(6) — arctan b), which can be found only by solving the joint system of Eq. (11) and that for rotational motion of the projectile around its mass center C.

4. Rotational Equation in Projective Coordinates

The equation (2) of the plane rotation around c. m. requires transition from time t to projective dual variable b according to the formulas

d0 d0 • . d20 • 2 d0 -w = d0 = d0t ■ = di0» +d50 (13>

The derivatives of the angle 0 by slope b are

d0_dtf 1 2b ~dh = ~db + 1 + b2, ~db'2 = db2 — (1 + b2)2. ( )

For the second derivative of slope b by t it receives after substitution of (10) in (7) and elementary transformations

.. g2 (ao(1+ e sin2 ^Vl + V ■ a%b + 2gyo sin 0(1 + b2) f) + d (70 sin 0(1 + b2)) f)

b =---T- .

a'bb u^ + 9^0 sin 0(1 + b2)3 \ bb /

Then it receives two ratios — the dimensionless one

„,2 -9(ao(l + e sin2 0)^1 + b2 ■ a'^J-1 + 2gl0 sin 0(1 + b2)3/2\ +

( b)2 \ \abb J

+ d(70 sin 0(1 + b2)3/2) (16)

second one

y(&)2 _ g(1 + b2)

. ___i_)

' V (abb)2 w- + no sin 0(1 + I,2)3/2

^ bb '

b2 + gio sin 0(1 + b2)3/2' 1 (a'')2 At + gjo sin 0(1 + b2)3/2

ub \ bb '

( a'b)2(1 + b2). (17)

The substitution (14)-(15) in the starting rotational equation (2) and dividing both its sides by b2 with taking into account (16)-(17) results in the next system

b b

= 2g (a0 (1 + e sin2 0(b)) + i0b sin 0(b)^j x

(a'bb)2

x a"bf Jr + 97o sin 0(1 + b2)3/2^ + b2 +g

1 , ■ ^Y^TÛ , Alo sin0(1 + b2)3/2)

b b d

d20 2b /d0 1 \ + Vd6 + 1 + b 2 )

d 2 (1 + 2)2 d 1 + 2

x a"^ g(ao{1 + £ sin2 0(b)) + lob sin 0(6)^^1 + b2 x (18)

x abb(\ + 2gio sin0(1 + b2)3/2) + ab b

+ io(1 + b2)3/2 cos0d0 + 3io sin 0 ■ 6(1 + 62)1/2+

+ T~ ■ (1 + b2) \ -T+--g9(1+ m3/2) - ~h(1 + b2)(a'lb)2 sin0.

JCz \ -r + gio sin 0(1 + b2)3/2 Jcz

It has total order of k = 3, i.e. two units less than nominal one of primal system (1)-(2) and one unit less of its real order. Besides, substituting r](b) = sin 0(6) we achieve its full algebraization up to square roots and exclude even such elementary functions as sin / cos 0. The last is more than important for numerical integration.

5. The Cauchy Problem in Projective-Dual Variables

After introducing into consideration the new dependent value

$(6) = + 910 sin 0(&)(1 + b2)3'2

abb(0)

which mechanical sense is the ratio —;r (5) the system (18) transforms to

x2

' d$(6)

db

-2g(a0 (l + £ sin2 §(b)) + y0b sin §(b)^j x

x f1+ gvo sin ^(&)(1 + b2)3/2 \ x/YTb2

V+$(&) - gyo sin #(&)(! + 62)3/V + '

d2$ 2b fd® 1

+ -;r +

- f— + d&2 (1 + b2)2 ' Vd6 + 1 + b2)

x ^g ^«o (l + e sin2 #(&)) + yob sin #(&)) \/l + &2 x $(&) + g7o sin tf(&)(1 + b2)3/2

(19)

+

+

($(6) - g7o sin 0(&)(1 + b2)3/2Y 7o(1 + b2)3/2 cos 0f + 37o sin 0 ■ 6(1 + b2)1/2 + ■ (1 + b2) ■ yQ^

a(1 + b2) sin

JCz ($(6) - gyo sin 0(&)(1 + b2)3/2)

2

Highlighted are the aerodynamic parameters to be determined experimentally in tube or somehow else.

As for initial conditions (ICs) for slope angle of the axis I and its time derivative it is assumed the possibility of some starting attack angle but not any rotation

6(&o) = Oo + #o,

di d& v o;

dQ(6o) db

dti(bo) _ d(Q - 0) _ 1 d& = d& = -1 + b2o'

For 6(6o) = 0 (6), then —= 0, hence

do

With account of (5) the initial launch condition of non-rotating projectile may be written as g

$(M = T/2 g 2Q , V(bo) = #o, <(M = - cos2 do. (20)

Vq cos2 vo

As for initially rotating projectile, say, the knife thrown rotating the initial conditions should take into account the equality (6) for 6(60)-value.

Thus to receive so-called resolventa-function [4] f (b) = abb and attack angle $(b) behavior on the trajectory it's necessary to solve the systems of n = 2 ODEs, one of second, another of the first order. Alternatively, after standard substitution Q(b) = $'b(b) the same as the systems of n = 3 eqs each of order k = 1. It follows that Legendre transformation decreases the dimensionality on one unit.

The solution of Cauchy problem (19)-(20) determines the coordinates and time as

x(&) = - i_^_

( ) J $(&0 - 9io sin0(60(1 + b'2)3/2

= 'I

bo

d

) - 37o sin0(6')(1 + 6/2)3/2;

) - gio sin0(6')(1 + 6/2)3/2'

(21)

6. The Wind Affection

Let it be tail- or headwind with the constant velocity u. It may be taken into consideration through simple recalculating of these two values: a) initial conditions for $(b) and b itself into reference frames connected with blowing wind in a standard way

« + a> sin do

po = tan vo =

Vo cos do — u

(22)

b) initial velocity with respect to the air as W0 = \/V0 + u2 — 2V0 cos 90u.

Then comeback transition is due to formulas

tan o =

Wo sin e'o

Wo cos d'o + u

Vo =

+ u2 + 2Wou cos 6'o.

(23)

After integrating the system (19) with thus recalculated ICs the absolute horizontal coordinate is recalculated as xabs = x(b) + ut(b) with x(b) and t(b) according to (21). As for vertical coordinate y(b) it is calculated by unchanged formula and the trajectory is determined parametrically. The absolute inclination is determined by this as

babs(b) =

y'(b)

x' (b) + ut '(b)'

As for perpendicular to initial velocity Vo side wind with speed w, its effect may be taken into account by recalculating in connected with the wind frames the initial inclination angle and start velocity Wo. The first is done according to following formula

^o = tan 6'0 =

V sin

cos2 e o + w2

tan 1

w

2 V2 cos2

w < Vo cos 0o, (24)

the second one is simply W^V2 + w2.

After the integration of (19) with the modified ICs the perpendicular coordinate is determined as zabs(b) = wt(b), the horizontal one is xabs(b) = \Jx2(b) — z2hs(b) ~ (b)

x(b) — W2x(b)' and the vertical yabs(b) = y(b) with x(b), y(b) and t(b) also according to (21).

As for arbitrary wind direction tp with respect to the axis Ox and hence the velocity of U = (u, 0, w) = (U cos p, 0, U sin p) the appropriate solution is received by applying twice the for-mulas above with account of small size of the perpendicular component w.

b

7. Numerical Integration, Main Results

It was executed by use of Maple 15 for the system of n = 3 ODEs of k = 1 order obtained by introduction of additional simple equation of $'b(b) = This results

in that one of two other nonlinear equations contains all three unknowns and another only two.

In such the system could be used easily to calculate the trajectory of sporting or combat projectile just after preliminary determination of aerodynamic coefficients and parameters.

Varied are the next values: initial velocity V0 = 0.15 + 0.7 and 1.5 + 2.3 of Mach, angle of throwing d0 = 0...85°, coefficients a0, e, 70, a and 5 for the forces and torques involved in wide range.

The next types of $(0)-behavior were found in different ranges of the parameters above.

1. A monotonous increase and the temporary stabilization of the fall velocity at about a flat maximum and for attack angle of up to about d0 + 2 (Fig. 2). All these take place for little ratios a/Jcz and different coefficients of damping 5. At abnormally large negative slopes stabilization is replaced by fall of attack angle

up to 90° what entails loss of speed along with increasing of the d . However to reveal this is possible only when throwing projectile from great height. As for the stabilization, it is not simple conservation of axis I direction in space for angular velocity w = djj yet differs slightly from zero and the attack angle is not anytime close to d0 + ^ which is expected for almost vertical landing. It depends in particular on damping coefficient S (thick and thin lines). However, at abnormally high its values the axis orientation conserves in fact. Actually this stabilization effect may occur, such as when shooting from a special arch or heavy crossbow arrows, poorly oriented to the velocity vector due from large moment of inertia.

2. Non-monotonous $(0)-dependence modulated by fast damping 6-oscillations of low "frequency" (Fig. 3) at mean —ratios and large 5, these oscillations being without damping in the limit S ^ 0 (Fig. 4).

3. High frequency damped 6-oscillations at large a/Jcz-ratios with decrement depending on 5 too (Fig. 5, 6).

V = 100, a0 := 0,0001, e = 10, 70 := 0,00005, Jz := 0,001, a := 3 x 10"9, S := 5 x 10"9 (thin black), 5 x 10"5 (thick black) 90 = 60°

V = 100, ao := 0,0001, e = 10, 70 := 0,00005, Jz := 0,001, a := 0,000003, S := 0,5 x 10"7 , 80 = 60°

\ /'

-12 -10 -2 0

-0,01- /

-0,02- /

-0,03c /

I-i?/tt_0,ld@/dt --- 0,0001F(ft) I

Figure 2. Attack angle 0(6), angular V(b) velocities vs b for small a/Jcz-ratios

^ and linear

at

Figure 3. This for mean a /Jcz and large 5

V = 100, a0 := 0,0001, e = 10, 70 := 0,00005, Jz := 0,001, a := 0,000003, 5 := 0, 0O = 60°

Figure 4. Attack angle $(b), angular w = and linear V(b) velocities vs

slope b at average a/Jcz-ratios and

S = 0

V = 100, a0 := 0,0001, e = 10, 70 := 0,00005, Jz := 0,001, a := 0,00003, S := 2 x 10"8 , 60 = 60°

-#/tt_0,01 de/dt--- 0,00001^(6)

Figure 5. This for large (r/Jcz and small 5

In general, beyond the vertex region of trajectory both attack angle and angular velocity behave in ¿-representation like classical oscillatory system with damping.

The 0(6) large varying, with no doubt, greatly affects the mass center trajectory, and this influence is defined mainly by the coefficient a of static torque. The less the a the wider range of 0(6) varying and the trajectory is shorter and lower (Fig. 7). Vice versa, the greater this coefficient the less amplitude of attack angle oscillations and the trajectory is closer to that of a heavy mass point. Other factors like lift power and damping coefficients in their real ranges are not able to affect the trajectory in such extent.

V = 100, ao := 0,0001, e = 10, 70 := 0,00005, Jz := 0,001, a := 0,00003, S := 2 x 10"7, 0O = 60°

0,001-

V = 100, a0 := 0,0001, e = 10, 70 := 0,00005, Jz := 0,001, a \= 5 x 10"8, 0O = 60°, a := 0,000003 (1), a := 0,0000003 (2), a ■- 5 x 10"9 (3), a := 5 x 10"8 (4)

-I?/tt

0,01 de/dt--- 0,00001F(6)

Figure 6. The ê(b) and df vs b at mean a/Jcz and large 5

Figure 7. The trajectory plots at different a/Jcz-ratios

It is worth to pay attention on essentially non-monotonous behavior of the speed at descending part of the trajectory when launching at high angles of throwing close to 90° (Fig. 8). This affect may be explained simply by projectile rotation arising at the top of trajectory from sudden aerodynamic impact due to reverse. And finally the attack angle is stabilized at 0 « k, 0 (mod 2^), e.i. after of about two turns.

e=5,V0=50,00=85 °, aO := 0.00050, Jz ■= 0.001, o := 0.000015, yO := 0.0002, S -5.10"7

V0=50mps,90=45deg,a0=0.0001,a=0.000005;Jz=0.001,E=10,Y= 0.00005

-120 -100 -80

-6-

x

|---a/ii-0.1 V(b) ---O.5d0/dt|

+10 mps -5 mps-0 — +5 mps

----lOmps_

Figure 8. Non-monotonic velocity behavior V(b) due to rotation when descending

Figure 9. Trajectory plots for the arrow launched at 60 = 45° with Vo = 50 mps by different tail- or headwind speeds

'0

As for the wind affection on the trajectory it is demonstrated by the plots below calculated for different wind velocities u = -10 ■ ■ ■ + 10 mps (Fig. 9) and the same projectile and launch parameters.

Thus, the developed method using dual-projective coordinates is more convenient than standard way with Cartesian coordinates and time. First it gives though para-metrically the trajectory equation and excludes non-important time variable.

Also it enables to fulfill numerically the qualitative analysis of such dynamic systems as systems with variable dissipation [3] and build up its phase portrait in w — 6 and $ — 6 planes for detailed study.

As for practical applications, the ODE system derived may use as an alternative method for verifying other ones in exterior ballistics of both sporting and combat projectiles. In addition, its application would allow improving these projectiles for better target engagement and more easily use. Finally the method when being elaborated may be implemented into ballistic calculators.

1. W. R. Rheingans, Exterior and interior ballistics of bows and arrows — review, Archery Review (1936) 236.

2. D. N. Gkritzapis, D. P. Margaris, E. E. Panagiotopoulos, et al., Prediction of the impact point for spin and fin stabilized projectiles, in: WSEAS Transactions on Information Science and Applications, Vol. 5, 2008, pp. 1667-1676.

3. M. V. Shamolin, Dynamical systems with variable dissipation: Approaches, methods, and applications, Fundamental and Applied Mathematics 14 (3) (2008) 3-327, in Russian.

4. V. V. Chistyakov, On one resolventa method for integrating the low angle trajectories of a heavy point projectile motion under quadratic air resistance, Computer Research and Modeling 3 (3) (2011) 161-171, in Russian.

8. Conclusions

References

5. V. V. Chistyakov, On integrating the projectile motion equations of a heavy point in medium with height decreasing density, Bulletin of Udmurt state university. Series "Mathematics. Mechanics. Computer Sciences" 1 (2012) 120-132, in Russian.

6. V. V. Chistyakov, Numerical-analytical integrating the equations of a point mass projectile motion at the velocities close to sonic peak of air drag exponent, Computer Research and Modeling 5 (5) (2013) 785-798, in Russian.

7. R. L. McCoy, Modern Exterior Ballistics: The Launch and Flight Dynamics of Symmetric Projectiles, Schiffer Publishing, Ltd., 1999.

УДК 531.55

Об одном методе численного интегрирования динамических уравнений плоскопараллельного полёта спортивного или боевого снаряда в условиях воздействия ветра

В. В. Чистяков

Инженерный факультет ФГБОУ ВПО «Ярославская государственная сельскохозяйственная академия» Тутаевское шоссе, 58, Ярославль, Россия, 150042

Стандартный путь интегрирования динамических уравнений для плоскопараллельного резистивного движения твердого тела подразумевает введение двух декартовых переменных x(t) и y(t) и угла атаки $(t) и, соответственно, трёх взаимосвязанных обыкновенных дифференциальных уравнений (ОДУ), каждое номинально II-го порядка. Это приводит к большому вычислительному объёму и рискам в точности получаемых решений. Предлагаемый метод исключает временную переменную t и уменьшает число функций до п = 2: угол атаки $(Ь) и подкасательная к траектории а(Ь), где b = tg в, а

в — угол наклона к горизонту вектора скорости V центра масс снаряда. Этот базирующийся на преобразованиях Лежандра подход делает интегрирование контролируемым и удобным особенно в рассматриваемом случае квадратичных по скорости аэродинамических усилий: лобовое сопротивление, подъёмная сила, консервативный и диссипатив-ный моменты. Также метод позволяет получить легко и надежно траектории снаряда в условиях встречного, попутного или бокового ветров. Исследованы основные области аэродинамических параметров, в которых имеет место различное поведение угла атаки $(Ь): квазистабилизация и апериодические автоколебания. Также обнаружено существенно немонотонное поведение величины скорости на участке падения с двумя минимумами при высоких углах запуска. Развитый метод может быть внедрён в процесс совершенствования реальных спортивных и боевых снарядов, таких как стрела лука, копьё, неуправляемый оперенный снаряд и др.

Ключевые слова: свободное резистивное движение, траектория, квадратичное сопротивление, подъёмная сила, консервативный и диссипативный моменты, угол атаки, проективно-двойственные переменные, ветер.

Литература

1. Rheingans W. R. Exterior and Interior Ballistics of Bows and Arrows — Review // Archery Review. — 1936. — P. 236.

2. Prediction of the Impact Point for Spin and Fin Stabilized Projectiles / D. N. Gkritzapis, D. P. Margaris, E. E. Panagiotopoulos et al. // WSEAS Transactions on Information Science and Applications. — Vol. 5, No 12. — 2008. — Pp. 1667-1676.

3. Шамолин В. А. Динамические системы с переменной диссипацией: подходы, методы, приложения // Фундаментальная и прикладная математика. — 2008. — Т. 14, № 3. — С. 3-327.

4. Чистяков В. В. Об одном резольвентном методе интегрирования уравнений свободного движения в среде с квадратичным сопротивлением // Компьютерные исследования и моделирование. — 2011. — Т. 3, № 3. — С. 265-277.

5. Чистяков В. В. Интегрирование уравнений свободного движения тяжёлой точки в среде с вертикальным градиентом плотности // Вестник Удмуртского университета. Серия «Математика. Механика. Компьютерные науки». — 2012. — Т. 1. — С. 120-132.

6. Чистяков В. В. Численно-аналитическое интегрирование уравнений свободного движения тяжелой точки вблизи звукового пика показателя степенного сопротивления // Компьютерные исследования и моделирование. — 2011. — Т. 3, № 3. — С. 785-798.

7. McCoy R. L. Modern Exterior Ballistics: The Launch and Flight Dynamics of Symmetric Projectiles. — Schiffer Publishing, Ltd., 1999.