An Integrated Program to Simulate the Multibody Dynamics of Flapping Flight

Nội dung text: An Integrated Program to Simulate the Multibody Dynamics of Flapping Flight

1. JST: Smart Systems and Devices Volume 31, Issue 1, May 2021, 076-083 An Integrated Program to Simulate the Multibody Dynamics of Flapping Flight Anh Tuan Nguyen 1*, Thanh-Dong Pham 1, Cong-Truong Dinh2, Jae-Hung Han3 1 Le Quy Don Technical University, Hanoi, Vietnam 2 Hanoi University of Science and Technology, Hanoi, Vietnam 3 Korea Advanced Institute of Science and Technology, Republic of Korea *Email: anhtuannguyen2410@gmail.com Abstract This paper presents an in-house developed program that couples multibody dynamics and aerodynamics codes to simulate flapping flight of insects and micro air vehicles. The multibody dynamics code is built based on the numerical solution of the Lagrange equation, while the extended unsteady vortex-lattice method is employed to develop the aerodynamics code. The solution from the governing equation is obtained by the use of the fourth-order Runge-Kutta method and validated against the simulation results from a commercial software MSC Adams for a micro air vehicle model. In this work, parallel computing techniques are applied while estimating the aerodynamics force to minimize the running time of the program. Keywords: Flapping flight, multibody dynamics, micro air vehicles, motion-aerodynamics interaction. 1. Introduction1 flying systems are simplified as only the six DOFs of the body are considered. Some authors have indicated The flapping flight of animals in the wild has that this simplification technique is valid only for been optimized by natural selection, which has lasted high-frequency flapping-wing systems with low wing for hundreds of millions of years (Ellington, 1991). mass. Due to these reasons, it is necessary to develop Nowadays, scientists have studied and adopted some an efficient approach that is capable of modeling the important features of insect and bird flight to improve multibody dynamics of flapping flight and the performance of flapping-wing micro-air vehicles considering the interaction between aerodynamics (FWMAVs) (Keennon and Klingebiel, 2012). In fact, and motion. flapping flight has some characteristics that do not exist in conventional fixed-wing and rotary-wing flights. While animals and micro-air vehicles flap xb their wings, the surrounding fluid moves in a such complex way that lift is generated to maintain a a) β balanced flight. Moreover, the dynamics of flapping flight is characterised by the multi-degree-of-freedom motion, including six degrees of freedom (DOFs) for xG the body, and three DOFs for each wing. The interaction of aerodynamics and motion in flapping zb flight causes lots of trouble to model it. Some zG Center of mass attempts have been made to simulate the multibody of the body dynamics of flapping flight while considering the two-way interaction between aerodynamics and motion. Zhang and Sun (2010) coupled the Navier- xwl xwr Stokes equations and the equations of motion for b) xb insect flight. However, this type of approach suffers from an excessive computation cost due to the use of y y a high-order aerodynamic model. Moreover, the wl wr yb solutions from this simulation program have not been xG rigorously validated. Another approach, which ignores the contribution of wing inertia, has been widely used to simulate and study flapping flight [2]. Using this approach, the dynamics properties of yG Fig. 1. The FWMAV model and the coordinate ISSN: 2734-9373 systems used in this study Received: 22 February 2020; accepted: 12 March 2021 76
2. JST: Smart Systems and Devices Volume 31, Issue 1, May 2021, 076-083 2.2. Aerodynamic Model The aerodynamic model is based on the extended UVLM [6], which was developed using the potential flow theory [3]. The present aerodynamic model has intermediate cost and fidelity and is capable of handling the nonlinearity of the flow. Some important effects of flapping-wing aerodynamics, such as the wake capture and the leading-edge vorstex can be modeled in this work. To reduce the running time of α the program while predicting the aerodynamic force, parallel computing techniques are applied. The validity of the present aerodynamic model for flapping-wing aerodynamics has been confirmed in previous studies [5]. φ θ 2.3. Multibody Dynamics Model Fig. 2. Euler angles to define the wing orientation. The present FWMAV consists of three objects, The program developed in this paper couples the including the body and the two wings with 12 DOFs in extended unsteady vortex-lattice method (UVLM), total. The dynamics of the FWMAV is modeled using which is regarded as an intermediate cost and fidelity the Lagrange method. The kinetic energy of the body aerodynamic method and the multibody dynamics Tb and each wing Tw are expressed as solver. The validity of the multibody dynamics solver 11GTT G b bb TMb= (VV b) bb+ (ΩΩ b) I bb is tested in a comparison with the solution from the 22 MSC Adams software. In this work, the program is (2) 11bTT b b bb applied to an insect-like FWMAV in hover. TMw= (VV w) ww+ (ΩΩ w) I ww 22 2. Methodlogy Here, M is mass, I is the tensor of moment of 2.1. FWMAV Model inertia, V and Ω are the velocity at the centre of mass and the angular velocity, respectively. The subscripts b The FWMAV model used in this study has the and w denote the body and a wing, whereas the same mass and geometric parameters as the hawkmoth superscripts G and b refer to the ground-fixed and Manduca sexta. The body and each wing weigh 1485.0 body-fixed coordinate system. and 46.87 mg, respectively; the wing length is 48.5 mm. Here, we use four coordinate systems, including the The velocity of the wing in the body-fixed b ground-fixed, the body-fixed and the two wing-fixed coordinate system Vw is coordinated systems as shown in Fig. 1. In this figure, VBVb= G +×Ω bb r + V b (3) β is an angle between the stroke plane of the wings and w G→ b b b wb// wb the body axis. where B → is the rotation matrix that transforms The wing orientation relative to the stroke plane Gb the velocity in the ground-fixed to the body-fixed is defined by three Euler angles φ , θ and α (Fig. 2). coordinate system; Vw/b and rw/b are the velocity and The sweep angle φ defines the up and down motions the position vector of the wing mass centre relative to of the wings, the elevation angle θ is used to control the body mass centre, respectively. the for- and backward motions, and α is the rotation angle of the wings about their feathering axes. The The Lagrange equation for the present FWMAV more detailed definitions of these angles can be found can be written in the following form [8]: in the literature [6]. The variations of φ , θ and α are dT∂∂ T −=Qη expressed by harmonic functions as follows: dt ∂∂ηη (4) η =xyz,,, ΦΘΨ , , φφ= − 0 cos( 2 πft) where T is the total kinetic energy; Q is the θθ= 0 cos( 4 πft) (1) generalized force derived from the virtual work done π ααπ= − 0 sin( 2 ft) by external forces, including gravitational and 2 aerodynamic forces; x, y and z are the displacements of the body centre of mass along the corresponding axes where φ , θ and α denote the amplitudes of the 0 0 0 in the ground-fixed coordinate system; Φ, Θ and Ψ are sweep, elevation and rotation angles, f is the flapping respectively the roll, pitch and yaw angles that are used frequency. In this study, it is assumed that angle β is to define the orientation of the body according to the 3- held constant. 2-1 sequence of rotations. 77
3. JST: Smart Systems and Devices Volume 31, Issue 1, May 2021, 076-083 The governing equation of the current multibody ∂Φ T = 1 Cη _ wwM dynamics problem is formulated as ∂η   M(Φ,t) Φ+ H( Φ,, Φ tt) =Q( Φ ,, ) (5) T ∂Φ T D= M 1 BrT b B η _w w G→ b wc/1 b × M, H, and Q matrices are given in detail as ∂η follows: T dB T CC++ C D+ D Gb→ b x__ b x wr x _ wl x__ wr x wl rB+ wc/1 b × dt CC++ C D+ D  y__ b y wr y _ wl y__ wr y wl T  b ∂Φ T d r CC++ C D+ D 1 T wc/ b × z__ b z wr z _ wl z__ wr z wl  EMη = BB→ ++Φ M =  _w w ∂η Gb 12 C+ C DD ++ D  dt ϕ__wr ϕwl ϕϕ__b wr ϕ _wl  T dB + ++ T b 1 Cθ__wr C θwl DDθθ__b wr D θ _wl Br G→ b wc/ b × dt + ++ Cψ_wr C ψ _wl DD ψψ __b wr D ψ _wl  ∂Φ T T b 1 ddBVG→ b b T wc/ b EE+ M VB+ x__ wr x wl w wc/ b G→ b ∂η dt dt EE+ y__ wr y wl η = xyz,, + EEz__ wr z wl H =  E+ E +− E( FF ++ F) and ϕ_b ϕ _wr ϕ _wl ϕϕ __b wr ϕ _wl E+ E +− E( FF ++ F) T θ_b θ _wr θ _wl θθ __b wr θ _wl ∂Φ C = M 2 BrT b B η _w w 1/wc b× G→ b E+ E +− E( FF ++ F) ψ_b ψ _wr ψ _wl ψψ __b wr ψ _wl ∂η G T ∂Φ T F 2 T bb Q = M Br r B+ T G w 1/wc b×× wc / b 1  ∂η BB→ M 1 Gb =  Dη _ w ∂Φ T 2 T b M, H and Q matrices are obtained from B11 IBw equation (4). Equation (5) is solved by the fourth- ∂η order Runge-Kutta method. T dB 1 b + [rBwc/ b] G→ b × The coefficients related to the body are dt expressed as T b ∂Φ d [rwc/ b ] EM= 2 BBT × ++Φ η _1w w Gb→ 1 ∂Φ T ∂η dt C = M 1  η _ bb T b dBGb→ ∂η Br1/[ wc b ] × dt T  ∂Φ = 2 ΒΒT b DIη _bb11 T T ∂η dB1 bb + [rrBwc/ b] [ wc /1 b ] T ×× T b dt dB bb  T ΒΒ 1 + ∂Φ d ( 11I b ) [rVwc// b] wc b b × 2  dt E = Φ d [rwc/ b ] T η _2b T × b +  TTB1 [ rBwc/1 b ] b ∂η dt × ∂∂ΦΦd [rwc/ b ] dt T b +MM22Φ ++BV× + w T 2w  1/wc b ∂∂ηηb dt η= ϕθψ d [rwc/ b ]  Here xyz,,,,, ; M and I respectively T b × Br B+ b 1/[ wc b ]× 1 dV dt T b wc/ b denote the mass and the moment of inertia. For the Br1/[ wc b ] × dt T dB  body, the moment of inertia is determined about its T bb 1 Br1/[ wc b] [ r wc / b ] ××dt center of mass. B1 is the matrix that converts the derivatives of the three Euler angles ϕθψ,, to the ddBBTT   11b ++bbΩ IBw 1/ Iw wb  angular velocity in the body-fixed coordinate system, dt dt  and B is a transform matrix. The superscripts b and G ∂∂ΦΦTTddIIbb   + 22TTww++ΦΩ b refer to the body-fixed and ground-fixed coordinate BB1 12  B 1wb /  ∂∂ηηdt dt T    systems, respectively. Φ1 is ( xyz,,) and Φ2 is b T b dB dΩ  BI 1 + T b wb/ T 1 w BI1 w  (ϕθψ,, ) . dt dt  The coefficients related to the wings are given η= ϕθψ,, in the following equations. [r]× denotes the skew- The generalized coordinate vector Φ is symmetric matrix of r, which is used to represent a T cross product in a matrix form. Φ=[xyzΦΘΨ] (6) The validity of the multibody dynamics solver is examined in a comparison with the solution from the 78
4. JST: Smart Systems and Devices Volume 31, Issue 1, May 2021, 076-083 MSC Adams software for the same FWMAV. Here, 3. Simulation results arbitrary external forces are applied to the two wings The present program is applied to the FWMAV and the body. Moreover, there is a phase difference in hover. The wing kinematic parameters are similar to between the motions of the left and right wings. Fig. 3 those of a real hawkmoth (Willmott and Ellington, shows the displacements and the orientations of the 1997) (Table 1). The initial velocities and the angular body from the present program and from MSC Adams. velocities are 0, and the initial pitch angle is 39.8 deg. MSC Adams is multibody dynamics software, Running the simulation with the present program which can provide a simulation framework for for the hovering FWMAV we can obtain the flight flapping-wing flight. However, the integration of trajectory. Next, constant lateral and vertically unsteady aerodynamics solvers into the software downward gust disturbances with the magnitude of framework is challenging because of its large 1.0 m/s is applied to the FWMAV, and the responses computational cost. Nguyen et al. [8] have made an of the FWMAV to these disturbances are investigated. attempt to do this kind of integration. However, it should be noted that in that study, the aerodynamics Table 1. Wing kinematic parameters solver is not entirely coupled with the iterations of the dynamics solver. Instead, the aerodynamic force is Parameters Values calculated for each time step, and this force is φ0 (deg) 60.0 assumed to be constant during the iterations of the θ (deg) 10.0 dynamics solver. This study presents the first attempt 0 to entirely couple the unsteady aerodynamics solver α0 (deg) 57.3 with a multibody dynamics code. f (Hz) 26.1 MSC Adams uses the GSTIFF solver for the β (deg) 54.8 multibody dynamics simulation. A predictor-corrector The variations of the trajectory parameters within method is applied to obtain the solution to the the time course of one wingbeat stroke cycle dynamics equations of the system. corresponding to the cases with and without the gust Fig. 5 and Fig. 6 show the comparisons for data disturbances are shown in Fig. 6 and 7. on the right and left wings, respectively. 20 2 40 15 30 0 10 20 5 10 -2 0 0 x-displacement (mm) x-displacement (mm) y-displacement MSC Adams (mm) z-displacement Present program -5 -4 -10 0 0.02 0.04 0.06 0.08 0.1 0 0.02 0.04 0.06 0.08 0.1 0 0.02 0.04 0.06 0.08 0.1 Time (s) Time (s) Time (s) 4 2 4 2 1 2 0 0 0 (rad) (rad) (rad) Φ Θ Ψ -2 -1 -2 -4 -2 -4 0 0.02 0.04 0.06 0.08 0.1 0 0.05 0.1 0 0.02 0.04 0.06 0.08 0.1 Time (s) Time (s) Time (s) Fig. 3. The validation against the MSC Adams software for the displacements and orientations of the body 79
5. JST: Smart Systems and Devices Volume 31, Issue 1, May 2021, 076-083 60 6 ) 1500 2 ADAMS In-house MBD 4 1000 40 2 500 20 0 0 0 -2 -500 -20 -4 -1000 Right wing marker velocities(m/s) -40 -6 -1500 Right wing marker accelerations(m/s Right wing marker displacements(mm) 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.02 0.04 0.06 0.08 0.1 0.12 Time (s) Time (s) Time (s) 4 10 300 6 200 4 100 2 0 0 -100 -2 Right wing angular velocities(rad/s) -200 -4 Right wing angular acclerations(rad/s) 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.02 0.04 0.06 0.08 0.1 0.12 Time (s) Time (s) Fig. 4. The validation against the MSC Adams software for the displacements of a marker on the right wing and the angular velocity and acceleration of right wing about three axes of the ground-fixed coordinate system 60 8 2000 ) ADAMS 2 In-house MBD 6 1500 40 4 1000 2 500 20 0 0 0 -2 -500 -4 -1000 -20 -6 -1500 Left wing marker velocities(m/s) -40 -8 Left wing marker accelerations(m/s -2000 Left wing marker displacements(mm) 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.02 0.04 0.06 0.08 0.1 0.12 Time (s) Time (s) Time (s) 4 10 400 8 6 300 4 200 2 100 0 0 -2 -100 -4 Left wing angular velocities(rad/s) -200 Left wing angular acclerations(rad/s) -6 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.02 0.04 0.06 0.08 0.1 0.12 Time (s) Time (s) Fig. 5. The validation against the MSC Adams software for the displacements of a marker on the left wing and the angular velocity and acceleration of left wing about three axes of the ground-fixed coordinate system 80
6. JST: Smart Systems and Devices Volume 31, Issue 1, May 2021, 076-083 Fig. 6. Displacements in the cases of lateral, vertical Fig. 7. Euler angles of the body in the cases of lateral, and no gust disturbances vertical and no gust disturbances a) YG XG YG ZG ZG XG b) YG XG YG ZG ZG XG Fig. 8. The movements of the body and the wings in the cases of the lateral gust (a) and the vertically downward gust (b) 81
7. JST: Smart Systems and Devices Volume 31, Issue 1, May 2021, 076-083 Pascal YG XG Z G Fig. 9. Wakes and pressure difference in the case of no gust (up), lateral gust (middle) and downward gust (down). The left subfigures show the simulation in the mid of the downstroke, the right ones are corresponding to the mid of the upstroke It is seen that when the lateral gust is applied, Moreover, strong wing-wake interaction can be the FWMAV is deflected largely about the roll axis visualized in this figure, which is the source of the and the motion along the y axis is divergent. challenge while handling the coupling between the Similarly, in the case of the vertically downward gust, unsteady aerodynamics solver and the multibody the pitch angle decreases and the FWMAV moves dynamics code. downwards. 4. Conclusions From these data, it is possible to state that in The paper has presented an in-house program hovering flight, the FWMAV is dynamically developed for the simulation of flapping flight. The unstable. The movements of the body and the wings motion-aerodynamics interaction problem is handled of the FWMAV in the cases of the gust disturbances by the integration of the multibody dynamics and are illustrated in Fig. 8. It is seen that the lateral gust aerodynamic codes. The multibody dynamics solver disturbance causes a bank turn motion. On the other was validated against the MSC Adams software. hand, the insect slightly pitches down in the case of When applying the program to a hawkmoth-like the downward vertical gust disturbance. This flapping-wing micro air vehicle model, the results response has demonstrated the dynamic instability of showed that this air vehicle is dynamically unstable the system. and easily deflected due to gust disturbances. Fig. 9 shows the wakes and nondimensional Acknowledgments pressure difference on the wings in the mid of the downstroke and upstroke phases. Here, data are This research is funded by Vietnam National nondimensionalized by the reference dynamic Foundation for Science and Technology pressure, in which we use the mean wing-tip speed. It Development (NAFOSTED) under grant number is seen that the wake moves with the gust direction. 107.01-2018.05. 82
8. JST: Smart Systems and Devices Volume 31, Issue 1, May 2021, 076-083 References stability of the hawkmoth Manduca sexta. Bioinspiration & Biomimetics, 12(1), p. 016007. [1] Ellington, C. P. (1991). Aerodynamics and the Origin of Insect Flight. Advances in Insect Physiology, 23, pp. 171–210. [6] Nguyen, A. T., Kim, J. K., Han, J. S., and Han, J. H. (2016b). Extended unsteady vortex-lattice method for insect flapping wings. Journal of Aircraft, 53(6), pp. [2] Lee, J.-S., Kim J.-K., and Han J.-H. (2015). Stroke 1709-1718. Plane Control for Longitudinal Stabilization of Hovering Flapping Wing Air Vehicles. Journal of Guidance, Control, and Dynamics, 38(4), pp. 800– [7] Willmott, A. P., and Ellington, C. P. (1997). The 805. mechanics of flight in the hawkmoth Manduca sexta. I. Kinematics of hovering and forward flight. Journal of experimental Biology, 200(21), pp. 2705-2722. [3] Katz, J., and Plotkin, A. (2001). Low-Speed Aerodynamics: From Wing Theory to Panel Methods. Cambridge University Press, New York. [8] Wittenburg, J., (2008). Dynamics of multibody systems. Springer, Germany. [4] Keennon, M., and Klingebiel, K. (2012). [9] Zhang, Y. L., and Sun, M. (2010). Dynamic flight Development of the nano hummingbird: A tailless stability of hovering model insects: Theory versus flapping wing micro air vehicle. In 50th AIAA simulation using equations of motion coupled with Aerospace Sciences Meeting Including the New Navier-Stokes equations. Acta Mechanica Sinica, Horizons Forum and Aerospace Exposition, p. 588. 26(4), pp. 509–520. [5] Nguyen, A. T., Han, J. S., and Han, J. H. (2016a). Effect of body aerodynamics on the dynamic flight 83