Flight control simulations of a butterfly-like flapping wing–body model by the immersed boundary–lattice Boltzmann method

Accepted Manuscript

Flight control simulations of a butterfly-like flapping wing–body model by the immersed boundary–lattice Boltzmann method Yuichi Nakatani, Kosuke Suzuki, Takaji Inamuro PII: DOI: Reference:

S0045-7930(16)30115-3 10.1016/j.compfluid.2016.04.027 CAF 3166

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

4 January 2016 8 April 2016 10 April 2016

Please cite this article as: Yuichi Nakatani, Kosuke Suzuki, Takaji Inamuro, Flight control simulations of a butterfly-like flapping wing–body model by the immersed boundary–lattice Boltzmann method, Computers and Fluids (2016), doi: 10.1016/j.compfluid.2016.04.027

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • Flight control methods for a butterfly-like flapping wing-body model

CR IP T

are proposed.

• The positions of the center of mass or the center of lift are changed.

• The center of mass is moved by an additional weight attached to the model.

AN US

• The center of lift is moved by a short interruption of the wing motion.

AC

CE

PT

ED

M

• Both methods can control the attitude of the model.

1

ACCEPTED MANUSCRIPT

CR IP T

Flight control simulations of a butterfly-like flapping wing–body model by the immersed boundary–lattice Boltzmann method Yuichi Nakatania , Kosuke Suzukib , Takaji Inamuroa,c,∗ a

AN US

Department of Aeronautics and Astronautics, Graduate School of Engineering, Kyoto University, Kyoto 615-8540, Japan b Institute of Engineering, Academic Assembly, Shinshu University, Nagano 380-8553, Japan c Advanced Research Institute of Fluid Science and Engineering, Graduate School of Engineering, Kyoto University, Kyoto 615-8540, Japan

Abstract

Flapping flight of insects and bio-inspired micro air vehicles requires not

M

only the generation of aerodynamic forces for supporting the weight against gravity and for driving forward, but also the flight control for maintaining

ED

the balance while flying. In this study, we investigate the flight control of a butterfly-like flapping wing–body model numerically by using the immersed

PT

boundary–lattice Boltzmann method. First, we simulate the control of a pitching motion by moving a weight along the body. It is found that the

CE

pitching angle can be controlled by determining the position of the weight with a P control scheme. Second, we simulate the control of a rolling motion

AC

from a large initial disturbance on the rolling angle by moving a weight perpendicular to the body. It is found that the rolling angle can be controlled ∗

Corresponding author. Tel.: +81 75 383 3770; fax: +81 75 383 3774. Email addresses: [email protected] (Kosuke Suzuki), [email protected] (Takaji Inamuro)

Preprint submitted to Computers & Fluids

May 2, 2016

ACCEPTED MANUSCRIPT

by determining the position of the weight with a PI control scheme. Finally, we simulate the control of a pitching motion by interrupting the wing motion

CR IP T

shortly at the top and/or the bottom dead centers of the wing motion. It is found that the pitching angle is controlled by determining the interval of the interruption with a PID control scheme. In addition, it is found that the present control works well even for a large disturbance.

Keywords: Flapping flight, Flight control, Fluid structure interaction,

AN US

Lattice Boltzmann, Immersed boundary 1. Introduction

In flapping flight of insects such as flies and butterflies and of bio-inspired micro air vehicles [1, 2], it is required not only to generate aerodynamic forces

M

for supporting the weight against gravity and for driving forward, but also to utilize the flight control for maintaining the balance while flying. While many

ED

efforts have been made over the last 50 years for investigating the generation of aerodynamic forces, the study of the flight control has a short history, only

PT

the past 15 years or so [3]. In such a short history, however, the flight control has become one of hot topics in the study of flapping flight through many

CE

pioneering studies. Taylor and Thomas [4] investigated the dynamic flight stability and the flight control of desert locusts (Schistocerca gregaria) in a

AC

steady forward flight. In this study, the steady forward flight is regarded as a fixed-point equilibrium, and small disturbances from the equilibrium are analyzed based on the linearized equations of motion of the insect’s body using time-averaged aerodynamic force and torque measured by experiments. The above simplification is called as the averaged-model theory [3]. Many 3

ACCEPTED MANUSCRIPT

researchers have applied the averaged-model theory to studies of the longitudinal stability and the lateral stability in hovering flight of several insects

CR IP T

[5, 6, 7, 8, 9]. In these studies, the hovering flight is regarded as a fixedpoint equilibrium. Sun and Wang [10] and Zhang and Sun [11] performed an

analysis on the stability and controllability and the stabilization control of a hovering model insect by applying the averaged-model theory using com-

putational data. More recently, Cheng et al. [12] studied the stabilization

AN US

controls of hovering hawk moths by combining the passive stability and controllability, the sensory systems and the feedback control responses. This study treats the hovering flight as an equilibrium and applies the averagedmodel theory using experimental data, too. A few pioneering works such as Cheng et al. [12] and Ristroph et al. [13, 14] have investigated the integrated

M

flight control systems including the sensory system. The more detailed and comprehensive review on the flight control in flapping flight can be found in

ED

the review article by Sun [3].

As seen above, most of studies on the flight control are based on the

PT

averaged-model theory and treat the hovering flight or the steady forward flight as a fixed-point equilibrium. However, a more important feature of

CE

flapping flight is in accelerated motions such as sudden starting, halting, and turning. In the accelerated motions, a take-off from a resting state is especially important. Since the accelerated motions are totally different

AC

from a fixed-point equilibrium flight, the averaged-model theory used for the hovering flight or the steady forward flight cannot be applied to them. Therefore, for further investigation on the flight control in the accelerated motions, it is desirable to analyze a nonlinear time-variant dynamic system

4

ACCEPTED MANUSCRIPT

including fluid motion, wing motion, body motion, and a control response. From the viewpoint of fluid dynamics, the nonlinear time-variant dynamic

CR IP T

system includes a highly complicated fluid–structure interaction. In the past few years, computational fluid dynamics become a potential candidate for analyzing the nonlinear time-variant dynamic system owing to recent advances in computational methods [15, 16, 17, 18]. Actually, a few researchers

[19, 20] have simulated the free flight of insects such as flies and butterflies.

AN US

More recently, Wu et al. [21] have simulated the free hovering flight with flight control of a fruit fly by using a PID control of wing kinematics. The

authors also have numerically investigated the free flight of a 2D symmetric flapping wing [22, 23], of a dragonfly-like flapping wing–body model [24], and of a butterfly-like flapping wing–body model [25] by using an immersed

M

boundary–lattice Boltzmann method (IB-LBM) [26]. Especially, Minami et al. [24] and Suzuki et al. [25] have studied not only the generation of aerody-

ED

namic forces but also the control of the free flight of the flapping wing–body models in the take-off and the transitional motion to a steady forward flight.

PT

In the present study, in order to advance the research by Suzuki et al. [25], we further discuss the flight control of the butterfly-like flapping wing–body

CE

model by using the IB-LBM [26]. The butterfly-like flapping wing–body model flaps downward for generating the lift force and backward for generating the thrust force like an actual butterfly. Therefore, the wing–body model

AC

utilizes aerodynamic forces parallel to the wingtip path, i.e. the drag-based thrust [27, 28]. Since insects uses the drag-based thrust in brisk maneuvers [27], the motion of the butterfly-like flapping wing–body model should be highly unstable. A flight control method which can stabilize the motion

5

ACCEPTED MANUSCRIPT

of the butterfly-like flapping wing–body model is expected to be a useful strategy for controlling arbitrary flapping air vehicles. Therefore, the main

CR IP T

objective of this study is not to clarify the flight control of an actual butterfly but to propose some control methods for flapping air vehicles.

In the above objective, we mainly focus on the control of pitching motion. Since the instability in pitching motion has been widely observed in

flapping flight by insects independently of species (flies [6, 19, 14], moths [12],

AN US

and bees [5] in hovering flight, and locusts [4] and butterflies [20] in forward flight), the control of pitching motion should be essential for arbitrary flapping air vehicles. While Suzuki et al. [25] investigated the control of pitching

motion by flexing the body of the wing–body model like an actual butterfly, in this study we propose two other control methods based on a more prim-

M

itive concept, i.e. the change of the relative position between the center of mass (COM) and the center of lift (COL). The COL is the point at which

ED

the net lift force acts, and the COL in the front/rear of the COM induces nose-up/down motion of aircrafts. Therefore, the relative position between

PT

the COM and the COL generally determines the rotational motion of aircrafts. Unlike fixed-wing aircrafts, the flapping wings dynamically changes

CE

the position of the COL and consequently the relative position between the COM and the COL. If the COM can be moved actively with the movement of the COL, the rotational motion of the flapping wing–body model may be

AC

able to be controlled. Therefore, in this study, we consider a control method in which the COM of the wing–body model is moved actively by using a moving weight attached to the body of the wing–body model. Hereafter we will call this control method the COM controller. On the other hand, the

6

ACCEPTED MANUSCRIPT

COL cannot be moved actively, since its position changes due to the wing motion. However, the time-averaged position of the COL can be changed by

CR IP T

modifying the wing motion. The simplest way to modify the wing motion is to interrupt the wing motion shortly. Therefore, in this study, we consider

another control method in which the wing motion is interrupted shortly at the top and/or the bottom dead centers of the flapping motion. Hereafter

we will call this control method the COL controller. Obviously, the COL

AN US

controller is mechanically much simpler than the COM controller.

The paper is organized as follows. In Section 2, we explain the butterflylike flapping wing–body model. In Section 3, we describe the governing equations and parameters of the system. The computational methods and conditions are presented in Section 4, and results and discussions are pre-

ED

2. Wing–body model

M

sented in Section 5. We finally conclude in Section 6.

We modify the butterfly-like flapping wing–body model originally pro-

PT

posed by Suzuki et al. [25] for investigating the flight control in flapping flight. In the followings, we first introduce the basic components and the

CE

basic wing motion of the original wing–body model. Although they are the same as those described by Suzuki et al. [25], we describe them again in order

AC

to define some symbols. Then, we explain the COM and COL controllers, which are the modifications from the original wing–body model. 2.1. Basic components The wing–body model is composed of the same two wings and a body. The wing is infinitely thin and has a square shape with side of length LW , 7

ACCEPTED MANUSCRIPT

CR IP T

Wing

Lw

Body

AN US

Connecting point

Figure 1: Illustration of (a) the butterfly-like flapping wing–body model [25] and (b) the coordinate system fixed to the body and the directions of pitching and rolling motions.

that is, the aspect ratio defined as the wing span squared divided by the wing area [29] is AR = 2. In addition, the flexibility and the mass of the wings are

M

neglected. The body is described as a straight infinitely thin rod with length Lb . We set Lb = LW for simplicity. The body has uniform (line) density ρb .

ED

Therefore, the COM of the body is at the middle point of the body, and the total mass of the body is Mb = ρb Lb . The two wings are connected to the

PT

body in a manner such that the middle point of the wing root is located at that of the body (see Fig. 1a).

CE

Here, we define the coordinate system fixed to the body (see Fig. 1b). Let the coordinate system fixed to the body be (X, Y, Z), the X-axis be parallel

AC

to the body, and its origin be located at the center of the body (COB). We define the forward direction relative to the body as the positive direction of the X-axis, the right direction as the positive direction of the Y -axis, and

the downward direction as the positive direction of the Z-axis.

8

ACCEPTED MANUSCRIPT

2.2. Basic wing motion The wing–body model flaps downward for generating the lift force and

CR IP T

backward for generating the thrust force like an actual butterfly. The wing motion is a combination of the attacking motion and the flapping motion,

described by the angle of attack α(t) and the flapping angle θ(t) at time t, respectively, as follows:

AN US

   2π αm 1 + cos t+γ , α(t) = 2 T   2π θ(t) = θm cos t , T

(1)

(2)

where αm is the maximum angle of attack, θm is the flapping amplitude corresponding to the half of the stroke amplitude, T is the period of flapping

M

motion, and γ is the phase shift. For more detail of the wing motion, see the reference [25].

ED

In this study, we set θm = 45◦ , αm = 90◦ , and γ = π/2. Fig. 2 shows the wing motion. It should be noted that with this set of parameters the wing– body model can generate an enough lift force for supporting actual insects’

PT

weight [25].

CE

2.3. COM and COL controllers While Suzuki et al. [25] investigated the control of a pitching motion by

AC

flexing the body of the wing–body model like an actual butterfly, the control method was somewhat complicated. In this study, we propose simpler and more primitive control methods for flapping air vehicles, i.e. the COM and the COL controllers described in Section 1. We apply the COM and the COL controllers to three kinds of the flight motion of the butterfly-like flapping 9

ACCEPTED MANUSCRIPT

t = 0.1T

t = 0.2T

t = 0.3T

t = 0.5T

t = 0.6T

t = 0.7T

t = 0.8T

t = 0.4T

CR IP T

t = 0.0T

AN US

t = 0.9T

Figure 2: Illustration of the basic wing motion for θm = 45◦ , αm = 90◦ , and γ = π/2.

wing–body model: (i) the pitching motion control by the COM controller; (ii) the rolling motion control by the COM controller; and (iii) the pitching motion control by the COL controller. The directions of pitching and rolling

M

motions are shown in Fig. 1(b). In the control (i), a point mass which can move along the body is attached to the body. By moving the point mass ac-

ED

tively, the COM of the model can be moved forward and backward. Since in general the relative position of the COM to the COL determines the direction

PT

of the rotation of aircrafts, it can be expected that the pitching motion of the wing–body model can be controlled by the forward-and-backward motion of

CE

the COM. In the control (ii), a point mass which can move perpendicular to the body is added to the model. In the similar way to the control (i), the

AC

rolling motion of the wing–body model is expected to be controlled by the right-and-left motion of the COM. In the control (iii), the wing motion is interrupted shortly at the top dead center (TDC) and/or the bottom dead center (BDC) of the flapping motion, i.e. at the time when θ(t) = θm and/or θ(t) = −θm . Although the COL is moved dynamically forward and back10

ACCEPTED MANUSCRIPT

ward due to the wing motion, the time-averaged position of the COL can be modified by the COL controller. Because of the same reason as the control

CR IP T

(i), it can be expected that the pitching motion of the wing–body model can be controlled by the pause of the wing motion. Specific algorithms of the controls (i), (ii), and (iii) are explained in Section 3. 3. Governing equations

AN US

In this section, we describe the governing equations of the system where

the fluid motion and the wing–body motion interact with each other. It should be noted that while the governing equations for the fluid motion are the same as those described by Suzuki et al. [25], the governing equations for

3.1. Fluid motion

M

the wing–body motion are modified with the COM and COL controllers.

ED

The fluid motion around the wings and the body is governed by the

∇ · u = 0,

(3)

Du 1 = − ∇p + ν∇2 u, Dt ρf

(4)

CE

PT

incompressible Navier–Stokes equations as follows:

where u is the fluid velocity, p is the pressure, ρf is the density of the fluid,

AC

and ν is the kinematic viscosity of the fluid. We consider the air at room temperature 20◦ C as the fluid, i.e. ρf = 1.205 [kg/m3 ] and ν = 1.512 × 10−5

[m2 /s]. It should be noted that the gravitational term is not appeared in Eq. (4). This is because the pressure p includes the gravitational potential. The no-slip condition should be satisfied on the surface of the wing–body 11

ACCEPTED MANUSCRIPT

model, i.e. the fluid velocity must be equal to the velocity of the wings and the body.

CR IP T

In this study, we take the mean wing tip speed defined by Utip = 4θm LW /T as the characteristic flow speed, and the wing length LW as the characteristic

length. The governing parameter of the above equations is the Reynolds number Re given by Re =

(5)

AN US

3.2. Wing–body motion

Utip LW . ν

In this section, the equations of the body motion with the COM controller for the controls (i) and (ii) are described at first. Next, the modifications of the wing motion in the control (iii) with the COL controller are described.

M

The coordinate system (X, Y, Z) fixed to the body and the coordinate system (x, y, z) fixed to the space are described in Fig. 3. The coordinate system

ED

(X, Y, Z) moves translationally and rotationally relative to the coordinate system (x, y, z). The gravitational acceleration G = 9.807 [m/s2 ] acts in the negative direction of the y-axis. The detailed definitions of the coordinate

PT

systems are found in the reference [25].

CE

3.2.1. Body motion with the COM controller (i) In the control (i), we attempt to control the pitching angle by actively

AC

moving a weight along the body. Since the wing motion is prescribed relative to the body as shown in Section 2.2, only the body motion with the COM controller is shown here. Let a weight with mass m be attached to the body with mass Mb and located at a point X = `. It should be noted that in the

present control the total mass of the model is M = Mb + m. The weight 12

CR IP T

ACCEPTED MANUSCRIPT

Figure 3: The configuration of the coordinate system (X, Y, Z) fixed to the body and the

AN US

coordinate system (x, y, z) fixed to the space.

moves only along the X-axis. The relative position ` is positive when the weight is located in the forward direction relative to the COB and negative when it is located in the backward direction relative to the COB, and it

M

changes in the range of −0.5Lb ≤ ` ≤ 0.5Lb . The pitching angle of the body θpitch is defined as the angle between the X-axis fixed to the body and

ED

the x-axis fixed to the space, and it is positive for the nose-up motion (see Fig. 4). We attempt to control the pitching angle θpitch by actively changing

PT

` as explained below.

Here, we assume that the body moves only in the x and y directions and

CE

rotates only in the pitching motion for simplicity, and we denote the vectors for the position, the velocity, and the force for only x and y components, and

AC

the torque for the z component. Let the position of the COB be (xb , yb ). The motion of the body–weight system is described by the three independent variables xb , yb , and θpitch . The Lagrangian for the body–weight system is

13

CR IP T

ACCEPTED MANUSCRIPT

Figure 4: The body of the wing–body model with the COM controller (i), i.e. control

AN US

of the pitching angle θpitch by moving a weight along the body actively, viewed from the z-axis.

given by:

M 2 m 2 ) (x˙ b + y˙ b2 ) + (`˙2 + `2 θ˙pitch 2 2 + m`˙ {x˙ b cos(θpitch ) + y˙ b sin(θpitch )}

M

L=

+ m`θ˙pitch {−x˙ b sin(θpitch ) + y˙ b cos(θpitch )}

ED

1 2 + I θ˙pitch − M yb G − m`G sin(θpitch ), 2

(6)

where the dot notation is used for time derivative, and I is the moment of

CE

PT

inertia of the body about the Y -axis as follows: I=

1 Mb L2b . 12

(7)

Let the aerodynamic force acting on the wings and the body be (fx , fy )

AC

and the aerodynamic torque acting on the wings and the body around the COB in the z direction be Tz . As mentioned in Section 2, the mass of the wings is neglected in the present wing–body model, and therefore the aerodynamic forces generated by the wings act on the body through the connecting point. Supposing that the relative position ` of the weight is 14

ACCEPTED MANUSCRIPT

determined as an input, the Lagrange equations for the motion of the body–

(8)

CR IP T

weight system can be obtained as follows:   d ∂L ∂L − = fx , dt ∂ x˙ b ∂xb   d ∂L ∂L − = fy , dt ∂ y˙ b ∂yb ! ∂L d ∂L − = Tz . dt ∂ θ˙pitch ∂θpitch

(9)

(10)

From the above equations, we can obtain the velocity (x˙ b , y˙ b ) of the COB

AN US

and the pitching angular velocity θ˙pitch . The kinematic equations, which determine the position (xb , yb ) of the COB and the orthogonal transformation matrix S for the rotational motion of (X, Y, Z) relative to (x, y, z), can be

found in the reference [25]. If m = 0, the above equations correspond to the

M

equations of the body motion without controller.

In the control (i), we consider the proportional (P) control scheme for

ED

determining the relative position ` as follows: `(t) = Kp (θ0 (t) − θpitch (t)),

(11)

PT

where θ0 is the desired pitching angle of the body and Kp is the proportional gain. In order to suppress the increase in the pitching angle θpitch , we set

CE

θ0 (t) = 0◦ . It should be noted that in this study the relative position ` is limited in the range of −0.5Lb ≤ ` ≤ 0.5Lb .

AC

The governing parameters of Eqs. (8)–(10) are the non-dimensional mass

NM and the Froude number F r defined as: M , ρf L3W Utip Fr = √ . LW G NM =

15

(12) (13)

ACCEPTED MANUSCRIPT

Center of the body

CR IP T

Weight

Figure 5: The body of the wing–body model with the COM controller (ii), i.e. control of

the back side of the wing–body model.

AN US

the rolling angle θroll by moving a weight perpendicular to the body actively, viewed from

In addition, there are two tuning parameters, i.e. m and Kp for the control. 3.2.2. Body motion with the COM controller (ii)

M

In the control (ii), we attempt to control the rolling angle by actively moving a weight perpendicular to the body. Let a weight with mass m move

ED

only along the Y -axis and be located at a point Y = `. It should be noted that in the present control the total mass of the model is M = Mb + m. The

PT

relative position ` is positive when the weight is located in the right direction relative to the COB and negative when it is located in the left direction

CE

relative to the COB, and it changes in the range of −0.5Lb ≤ ` ≤ 0.5Lb . The rolling angle of the body θroll is defined as the angle between the Y -axis

AC

fixed to the body and the z-axis fixed to the space, and it is positive when the Y -axis is inclined right downward (see Fig. 5). We attempt to control the rolling angle θroll by actively changing ` as explained below. Here, we assume that the body rotates only in the rolling motion for

simplicity, while moving in the x, y, and z directions. Let the position of the 16

ACCEPTED MANUSCRIPT

COB be (xb , yb , zb ). The motion of the body–weight system is described by the four independent variables xb , yb , zb , and θroll . The Lagrangian for the

L=

CR IP T

body–weight system is given by: m M 2 2 ) (x˙ b + y˙ b2 + z˙b2 ) + (`˙2 + `2 θ˙roll 2 2 + m`˙ {−y˙ b sin(θroll ) + z˙b cos(θroll )}

− m`θ˙roll {y˙ b cos(θroll ) + z˙b sin(θroll )}

(14)

AN US

1 2 + Ir θ˙roll − M yb G + m`G sin(θroll ), 2

where Ir is the moment of inertia of the body about the X-axis. Since the body of the wing–body model is an infinitely thin rod, Ir should be zero. In this study, however, we assume that the wing–body model has a virtual

M

moment of inertia about the X-axis calculated as that of a circular cylinder with length Lb and diameter d. We set d = 0.16Lb , which is close to the

ED

data of a small butterfly named Janatella leucodesma [30]. Consequently, Ir

PT

can be calculated as follows:

1 Ir = Mb (0.16Lb )2 . 8

(15)

Let the aerodynamic force acting on the wings and the body be (fx , fy , fz )

CE

and the aerodynamic torque acting on the wings and the body around the COB in the x direction be Tx . Supposing that the relative position ` of the

AC

weight is determined as an input, the Lagrange equations for the motion of

17

ACCEPTED MANUSCRIPT

(16)

CR IP T

the body–weight system can be obtained as follows:   d ∂L ∂L − = fx , dt ∂ x˙ b ∂xb   ∂L d ∂L − = fy , dt ∂ y˙ b ∂yb   d ∂L ∂L − = fz , dt ∂ z˙b ∂zb   ∂L d ∂L = Tx . − ˙ dt ∂ θroll ∂θroll

(17) (18) (19)

AN US

From the above equations, we can obtain the velocity (x˙ b , y˙ b , z˙b ) of the COB and the rolling angular velocity θ˙roll . The kinematic equations can be found in the reference [25]. If m = 0, the above equations correspond to the equations of the body motion without controller.

In the control (ii), we consider the proportional–integral (PI) control

M

scheme for determining the relative position ` as follows: Z t `(t) = Kp (θ0 (t) − θroll (t)) + Ki (θ0 (t0 ) − θroll (t0 ))dt0 ,

(20)

ED

0

where θ0 is the desired rolling angle of the body and Kp and Ki are the

PT

proportional and integral gains, respectively. In order to suppress the increase in the rolling angle θroll , we set θ0 (t) = 0◦ . It should be noted that in this

CE

study the relative position ` is limited in the range of −0.5Lb ≤ ` ≤ 0.5Lb . The governing parameters of Eqs. (16)–(19) are the non-dimensional mass

AC

NM and the Froude number F r defined by Eqs. (12) and (13). In addition, there are three tuning parameters, i.e. m, Kp , and Ki for the control. 3.2.3. Body motion and wing motion with the COL controller (iii) In the control (iii), we attempt to control the pitching angle by interrupting the wing motion shortly at the TDC and/or the BDC. The equations of 18

ACCEPTED MANUSCRIPT

the body motion are not modified by the COL controller from those without controller. Therefore, the equations of the body motion with the COL

1 2 M 2 − M yb G, (x˙ b + y˙ b2 ) + I θ˙pitch 2 2

d dt

  d ∂L ∂L − = fx , dt ∂ x˙ b ∂xb   d ∂L ∂L − = fy , dt ∂ y˙ b ∂yb ! ∂L ∂L − = Tz . ∂θpitch ∂ θ˙pitch

AN US

L=

CR IP T

controller correspond to Eqs. (6)–(10) with m = 0, that is, (21)

(22)

(23) (24)

The total mass M is equal to the mass of the body Mb in this control. On

M

the other hand, the wing motion is modified by the COL controller from the basic wing motion shown as Eqs. (1) and (2). Let the wing motion be

ED

interrupted by a short interval Ttop at the TDC and/or by Tbottom at the BDC of the flapping motion. The intervals Ttop and Tbottom are determined as inputs. Therefore, the angle of attack α and the flapping angle θ are

AC

CE

PT

modified as follows:   αm  , (0 ≤ t < Ttop ),   2    
  αm 2π   1 + cos (t − T ) + γ , (Ttop ≤ t < Ttop + T2 ), top  2 T   α(t) = αm , (Ttop + T2 ≤ t < Ttop + Tbottom + T2 ), 2     
  αm  1 + cos 2π (t − Ttop − Tbottom ) + γ ,  2 T       (Ttop + Tbottom + T2 ≤ t < Ttop + Tbottom + T ), 19

(25)

ACCEPTED MANUSCRIPT

(26)

CR IP T

   θc + θm , (0 ≤ t < Ttop ),     
   θc + θm cos 2π (t − Ttop ) , (Ttop ≤ t < Ttop + T2 ),  T   θ(t) = θc − θm , (Ttop + T ≤ t < Ttop + Tbottom + T ), 2 2   
   θc + θm cos 2π (t − Ttop − Tbottom ) ,  T       (Ttop + Tbottom + T2 ≤ t < Ttop + Tbottom + T ),

where θc is the center of the flapping angle. Eqs. (25) and (26) are illustrated

AN US

in Fig. 6, and it can be found that both α and θ are continuous. If θc = 0 and

Ttop = Tbottom = 0, the above equations correspond to Eqs. (1) and (2) of the wing motion without controller. In the control (iii), since the wing motion is interrupted, the generated lift force should be reduced. In order to prevent the wing–body model from falling and colliding the ground, we enhance the

M

lift force by changing the center of the flapping angle θc when the vertical velocity of the body is quite small. If the center of the flapping angle θc is

ED

negative, i.e. the flapping motion is oriented downward, it is expected that the lift force is enhanced due to the downward-oriented flapping motion.

PT

In the control (iii), we consider the proportional–integral–derivative (PID)

AC

CE

control scheme for determining the input intervals Ttop and Tbottom as follows: Z t t t (θ0 (t0 ) − θpitch (t0 ))dt0 + Kdt (θ˙0 (t) − θ˙pitch (t)), Ttop (t) = Kp (θ0 (t) − θpitch (t)) + Ki 0

(27)

Tbottom (t) = Kpb (θ0 (t) − θpitch (t)) + Kib

Z

0

t

(θ0 (t0 ) − θpitch (t0 ))dt0 + Kdb (θ˙0 (t) − θ˙pitch (t)), (28)

where θ0 is the desired pitching angle of the body and Kpt , Kit , and Kdt are the proportional, integral, and derivative gains for Ttop , respectively, and Kpb , 20

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 6: Illustration of the wing motion described by Eqs. (25) and (26) in the control (iii).

Kib , and Kdb are the gains for Tbottom . In order to suppress the increase in

M

the pitching angle θpitch , we set θ0 (t) = 0◦ . It should be noted that the input interval Ttop and Tbottom are calculated by Eqs. (27) and (28) only at the

Tbottom ≥ 0 .

ED

TDC and the BDC, respectively, and limited in the ranges of Ttop ≥ 0 and The governing parameters of this system are the same as those without

PT

controller, i.e. the non-dimensional mass NM and the Froude number F r defined by Eqs. (12) and (13). In addition, there are seven tuning parameters,

CE

i.e. Kpt , Kit , Kdt , Kpb , Kib , Kdb , and θc for the control.

AC

3.3. Governing parameters As shown in Sections 3.1 and 3.2, the governing parameters of the system

are the Reynolds number Re, the non-dimensional mass NM , and the Froude number F r for the three controls (i), (ii), and (iii). It should be noted that the tuning parameters for the controls such as m, Kp , Ki , and Kd are not 21

ACCEPTED MANUSCRIPT

included in the governing parameters. In order to calculate free flights of the wing–body model in fluid, we have to determine the three parameters. In this

CR IP T

study, we set the governing parameters and the corresponding dimensional parameters as shown in Table 1. This set of parameters is close to that of a

fruit fly, which is a typical flapping flyer. If the control methods presented above are available for this set of parameters, it is expected that they can be applied to flight control for artificial flapping air vehicles such as insect-scale

AN US

MAVs.

Table 1: The governing parameters and the corresponding dimensional parameters used in this study.

F r/Re

NM

300

2.9 × 10−2

38

ρf (kg/m3 )

ν (m2 /s)

1.205

1.512 × 10−5

M

Re

G (m/s2 ) LW (mm) Utip (m/s) M (mg) 3.0

1.512

1.236

ED

9.807

PT

3.4. Degrees of freedom

CE

In the control (i), we assume that the body moves only in the x and y directions and rotates only in the pitching motion. Hence, the motion of the wing–body model has 3 degrees of freedom (DOFs). In the control (ii), we

AC

assume that the body rotates only in the rolling motion, while moving in the x, y, and z directions. Thus, the motion of the wing–body model has 4

DOFs. In the control (iii), we assume that the body moves only in the x and

y directions and rotates only in the pitching motion. Therefore, the motion

22

ACCEPTED MANUSCRIPT

of the wing–body model has 3 DOFs. Controlling all rotational motions with

4. Computational methods and conditions

CR IP T

6 DOFs is remained in future work.

In this study, we use the IB-LBM [26] and the second-order Adams–

Bashforth method for calculating the fluid motion and the wing–body motion, respectively, in the same way as Suzuki et al. [25]. The fluid motion,

AN US

the wing–body motion, and the control responses such as `, Ttop , and Tbottom are calculated with a weak-coupling algorithm as follows:

1. The motion (translational and rotational velocities) and the configuration (position and attitude) of the wing–body model are updated by

M

one time step by using the aerodynamic forces and the control responses calculated in the previous time steps.

ED

2. The fluid velocity and the pressure are updated by one time step so that the no-slip boundary condition is satisfied on the surface of the wing–body model with the updated motion and configuration.

PT

3. The aerodynamic forces acting on the wing–body model are calculated from the updated flow flied.

CE

4. The control responses are determined by using the updated attitude of the wing–body model.

AC

5. One time step is advanced.

For the detail of the numerical method, see the reference [26]. The computational conditions are the same as those described in the

reference [25]. The computational domain is a cube with width W = 12LW 23

ACCEPTED MANUSCRIPT

No slip

CR IP T

L

y

L L

x

AN US

z

Figure 7: Computational domain for simulations of flows around the butterfly-like flapping wing–body model.

as shown in Fig. 7. The boundary condition on two sides perpendicular to the

M

x-axis is the periodic boundary condition, and on the other sides the no-slip condition is used. The inner fine grid [31] used only around the wing–body

ED

model is a cube with width D = 2.4LW . The COB is initially placed at the center of the domain filled with a stationary fluid at uniform pressure. In

PT

this study, we set LW = 50∆x and T = 6000∆t, which is the same condition as that used by Suzuki et al. [25] for Re = 300.

CE

5. Results and discussion

AC

In this section, we show some results of preliminary calculations as well

as the results for the controls (i), (ii), and (iii). All cases considered in the present study are listed in Appendix A.

24

ACCEPTED MANUSCRIPT

5.1. Pitching motion without control Before the flight controls with COM and COL controllers, we show the

CR IP T

result without control, which corresponds to those of the control (i) with m = 0 (see Section 3.2.1). Although it was presented in our previous study [25], we will show it again for helping readers to understand the necessity of flight control.

Fig. 8 shows the time variation of the pitching angle of the body, the

(b) 1.5

100 80

1

y b / LW

[deg]

(a)

AN US

trajectory of the COB in the x-y plane, and the attitude of the wing–body

60 40 20

0.5

0 0

2

4

6

t/T

10

-0.5 -0.5

0

0.5

1

1.5

2

2.5

3

x b / LW

CE

PT

ED

(c)

8

M

0

Figure 8: Results without control: (a) time variation of the pitching angle of the body,

AC

(b) trajectory of the COB, and (c) attitude of the wing–body model viewed from the right side of the model at t/T = 4, 6, and 8. In (b), the initial position of the COB is denoted

by (xb , yb ) = (0, 0), and the symbols on the trajectories indicate the position of the COB when the wings are at TDC. In (c), the wing–body model is shown as a red surface and the isosurface of the magnitude of the vorticity (|∇ × u| = 2Utip /LW ) is shown as a purple surface.

25

ACCEPTED MANUSCRIPT

model. It can be seen from Fig. 8(a) that the pitching angle θpitch increases monotonically. We can see from Figs. 8(b) and (c) that while the wing–

CR IP T

body model goes upward and forward, the model gradually raises its head and consequently gets off-balance. This means that the wing–body model is longitudinally unstable. Therefore, the control of the pitching motion should be essential. 5.2. Pitching motion with the COM controller (i)

AN US

First, we simulate the pitching motion with the COM controller (i), i.e. by actively moving a weight along the body of the wing–body model.

As a preliminary calculation, we simulate the free flight of the wing–body model with the weight of m = 0.2M fixed at the COB (` = 0), the tail of the

M

body (` = −0.5Lb ), and the head of the body (` = 0.5Lb ). Fig. 9 shows the time variation of the pitching angle of the body, the trajectory of the COB

ED

in the x-y plane, and the attitude of the wing–body model. It can be seen from Fig. 9(a) that the pitching angle θpitch increases monotonically for ` = 0 and −0.5Lb while θpitch decreases for ` = 0.5Lb . From Figs. 9(b), (c1), and

PT

(c2), we can see that for ` = 0 the wing–body model goes upward and raises its head in the same way as the results without control, while for ` = 0.5Lb

CE

the wing–body model goes downward and lowers its head. It suggests that the time-averaged position of the COL is located between ` = 0 and 0.5Lb ,

AC

and whether the pitching angle increases or decreases is determined by the position of the weight. It can be expected from this result that the pitching angle can be controlled by actively moving the weight. Next, we calculate the motion of the wing–body model when the relative position ` of the weight is determined by Eq. (11). In this simulation, we 26

ACCEPTED MANUSCRIPT

set m = 0.2M and Kp = −0.1LW /1◦ = −5.73 × LW as an example. Since

(a)

(b)

90

1.5

0

0.5 0 -0.5

-45

= 0.0 Lb = 0.5 Lb = - 0.5 Lb

-1 -90 0

2

4

6

8

-1.5 -0.5

10

t/T

CR IP T

y b / LW

[deg]

1 45

0

0.5

1

1.5

2

2.5

3

3.5

4

x b / LW

AN US

(c1)

PT

ED

M

(c2)

Figure 9: Results when the additional weight is fixed at the tail of the body (` = −0.5Lb ), the COB (` = 0) and the head of the body (` = 0.5Lb ): (a) time variation of the pitching

CE

angle of the body, (b) trajectory of the COB, and (c) attitude of the wing–body model with the weight viewed from the right side of the model for (c1) ` = 0 and (c2) ` = 0.5Lb at

AC

t/T = 4, 6, and 8. In (b), the initial position of the COB is denoted by (xb , yb ) = (0, 0), and the symbols on the trajectories indicate the position of the COB when the wings are at TDC. In (c1) and (c2), the wing–body model is shown as a red surface, the weight is shown as a green point, and the isosurface of the magnitude of the vorticity (|∇×u| = 2Utip /LW ) is shown as a purple surface.

27

ACCEPTED MANUSCRIPT

the above preliminary calculations suggest that ` should be greater than 0 for θpitch ≥ 0, we use negative Kp . Fig. 10 shows the time variation of the

CR IP T

pitching angle of the body, the trajectory of the COB in the x-y plane, the time variation of the relative position ` of the weight, and the attitude of

the wing–body model. Note that ‘Without Control’ in Figs. 10(a) and (b)

represents the result for m = 0, i.e. the same result as in Fig. 8. It can be seen from Figs. 10(a) and (b) that with the present control the pitching angle is

AN US

controlled in the range of |θpitch | < 10◦ up to 20 strokes, and the wing–body model can go upward against gravity. In addition, we can see from Fig. 10(c)

that the relative position ` of the weight oscillates with large amplitude in the early stage of the take-off, and then oscillates around ` = 0.2Lb with the amplitude about 0.2Lb . Due to the motion of the weight, the wing–body

M

model can go upward and forward stably as shown in Fig. 10(d).

ED

5.3. Rolling motion with the COM controller (ii) Second, we simulate the rolling motion with the COM controller (ii),

model.

PT

i.e. by actively moving a weight perpendicular to the body of the wing–body

We preliminarily simulate the free flight of the wing–body model with the

CE

weight of m = 0.2M fixed at the COB (` = 0), the left hand side of the body (` = −0.5Lb ) and the right hand side of the body (` = 0.5Lb ). Fig. 11 shows

AC

the time variation of the rolling angle of the body and the trajectory of the COB in the z-y plane. It can be seen from Fig. 11(a) that the rolling angle

θroll does not change for ` = 0 while it increases and decreases monotonically for ` = 0.5Lb and ` = −0.5Lb , respectively. Since the wing motion of the wing–body model is bilaterally symmetrical, the result for ` = 0 maintains 28

(a)

80 60

yb / LW

[deg] pitch

(b)

Without Control Controlled

100

40 20 0 0

5

10

15

4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 -0.5 -2

20

0

2

(c) 0.6

/ Lb

-0.2 -0.4 -0.6 0

6

8

10

12

AN US

0.4

0

4

xb / LW

t/T

0.2

CR IP T

ACCEPTED MANUSCRIPT

5

10

15

20

t/T

PT

ED

M

(d)

Figure 10: Results for the pitching motion control by the COM controller: (a) time

CE

variation of the pitching angle of the body, (b) trajectory of the COB, (c) time variation of the relative position of the weight, and (d) attitude of the wing–body model viewed from the right side of the model at t/T = 4, 6, and 8. In (b), initial position of the COB

AC

is denoted by (xb , yb ) = (0, 0), and the symbols on the trajectories indicate the position

of the COB when the wings are at the TDC. In (c), dashed lines denote ends of the body. In (d), the wing–body model is shown as a red surface, the weight is shown as a green point, and the isosurface of the magnitude of the vorticity (|∇ × u| = 2Utip /LW ) is shown as a purple surface.

29

ACCEPTED MANUSCRIPT

the bilateral symmetry, i.e. θroll = 0◦ . Next, we calculate the motion of the wing–body model when the relative

CR IP T

position ` of the weight is determined by Eq. (20). In this simulation, we attempt to maintain the balance of the wing–body model from an initial disturbance of the rolling angle. We set θroll (0) = −30◦ as an initial distur-

bance. As for tuning parameters for the control, we set m = 0.2M , Kp = −0.02LW /1◦ = −1.15×LW , and Ki = −0.012LW /(1◦ T ) = −0.688×LW /T as

AN US

an example. Fig. 12 shows the time variation of the rolling angle of the body, the trajectory of the COB in the z-y plane, the time variation of the relative position ` of the weight, and the vorticity around the wing–body model.

Note that ‘Without Control’ in Figs. 12(a) and (b) represents the result for

150

θ roll [deg]

100

0

-150

PT

-50

2

CE

0

4

6

1.5 1 0.5

ED

50

-100

(b)

= 0.0 Lb = - 0.5 Lb = 0.5 Lb

yb / Lw

(a)

M

m = 0, and the horizontal axis in Fig. 12(b) is not the forward displacement

0 -0.5 -1 -1.5 -2

8

10

12

14

-5 -4 -3 -2 -1

t/T

0

1

2

3

4

5

zb / Lw

AC

Figure 11: Results when the additional weight is fixed at the COB (` = 0), the left hand side of the body (` = −0.5Lb ), and the right hand side of the body (` = 0.5Lb ): (a) time

variation of the rolling angle of the body and (b) trajectory of the COB. In (b), the initial position of the COB is denoted by (zb , yb ) = (0, 0), and the symbols on the trajectories indicate the position of the COB when the wings are at the TDC.

30

ACCEPTED MANUSCRIPT

40

(b) 5

20

4 yb / LW

roll[deg]

0 -20 -40

Without Control Controlled

-60

1

-1

-100 0

5

10

15

20

25

30

-3 -2.5 -2 -1.5 -1 -0.5 zb / L W

0

0.5

1

AN US

t/T

(d)

0.6 0.4 0.2

/ Lb

2

0

-80

(c)

3

CR IP T

(a)

0 -0.2

-0.6 0

5

10

15

20

25

ED

t/T

M

-0.4

PT

Figure 12: Results for the rolling motion control by the COM controller: (a) time variation of the rolling angle of the body, (b) trajectory of the COB, (c) time variation of the relative

CE

position of the weight, and (d) vorticity around the wing–body model at t = 23.2T viewed from the front side of the wing–body model. In (b), initial position of the COB is denoted by (zb , yb ) = (0, 0), and the symbols on the trajectories indicate the position of the COB

AC

when the wings are at the TDC. In (c), dashed lines denote ` = −0.5Lb and 0.5Lb . In (d),

the wing–body model is shown as a red surface, the moving weight is shown as a green point, and the isosurface of the magnitude of the vorticity (|∇ × u| = 2Utip /LW ) is shown as a purple surface.

31

ACCEPTED MANUSCRIPT

but the lateral displacement of the COB. It can be seen from Fig. 12(a) that with the present control the rolling angle can be controlled in the range of

CR IP T

|θroll | < 5◦ after a large rolling motion in the early stage of the take-off, while without the control θroll increases monotonically. Fig. 12(b) shows that the

wing–body model with the control can go upward with a small lateral mo-

tion. In addition, we can see from Fig. 12(c) that the relative position ` of the weight oscillates with the large amplitude at first, and then it oscillates

AN US

around ` = −0.2Lb with the small amplitude. One might consider that the

weight should oscillate around ` = 0 when the steady flight. However, even at t/T = 23, the flow field around the model is laterally asymmetrical due to the history effect of flow (see Fig. 12d), and consequently the roll moment

cancel out the roll moment.

M

is induced. Therefore, the weight oscillates around ` = −0.2Lb in order to

ED

5.4. Pitching motion with the COL controller (iii) Finally, we simulate the pitching motion with the COL controller (iii), i.e. by interrupting the wing motion shortly at the TDC and/or the BDC.

PT

In this control simulation, we perform three preliminary calculations at first. Then, we conduct a free flight simulation with the COL controller. Finally,

CE

we simulate the response to a sudden large disturbance. As a first preliminary calculation, we calculate the time variation of

AC

the pitching torque when the equations of the body motion (22)–(24) are neglected (i.e. the body is fixed) for Ttop = Tbottom = 0 and θc = 0◦ . Fig. 13 shows the time variation of the pitching torque coefficient CM = 2 Tz /(ρf Utip L3W ) in one stroke of 9.0 ≤ t/T ≤ 10.0. We can see from Fig. 13

that the curve of CM has a shape similar to a square wave which is negative 32

ACCEPTED MANUSCRIPT

1

0

-0.5

-1

CR IP T

0.5

9 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10

t/T

AN US

Figure 13: Time variation of the pitching torque coefficient CM when the body of the wing–body model is fixed without controller.

in 9.2 ≤ t/T ≤ 9.7 and positive in the other duration. Since at the TDC (t/T = 9.0) and the BDC (t/T = 9.5) the pitching torque is positive and

M

negative, respectively, it can be expected that interrupting the wing motion at the TDC and the BDC induces the nose-up and the nose-down motion,

ED

respectively, due to the rotational fluid inertia. As a second preliminary calculation, we examine the effect of the center of

PT

the flapping motion θc on the translational motion of the body. We calculate the equations of the translational motion of the body (22) and (23) and

CE

neglect that of the rotational motion (24) for Ttop = Tbottom = 0. Fig. 14 shows the trajectory of the COB for θc = 0◦ , 20◦ and −20◦ . It can be

AC

seen from this figure that the wing–body model can go upward faster as θc decreases. This result means that the lift force can be enhanced by setting θc negative. As a third preliminary calculation, we simulate free flight of the wing–

body model when the wing motion is interrupted by a constant short interval

33

ACCEPTED MANUSCRIPT

6 5

3 2

CR IP T

yb / LW

4

θc = 0◦

1

= 20◦ = − 20◦

0 -2

0

2

4

6

8

10 12 14 16 18 20

AN US

xb / L W

Figure 14: Trajectory of the COB for the center of the flapping angle of θc = 0◦ and ±20◦ when the equation of the rotational motion is neglected for Ttop = Tbottom = 0. The initial position of the COB is denoted by (xb , yb ) = (0, 0), and the symbols on the trajectories

M

indicate the position of the COB when the wings are at the TDC.

Ttop = 0.2T at the TDC or by Tbottom = 0.2T at the BDC every flapping

ED

period. Fig. 15 shows the time variation of the pitching angle of the body and the trajectory of the COB in the x-y plane. It can be seen that the

PT

pitching angle increases when the wing motion is interrupted at the TDC, while it decreases when the wing motion is interrupted at the BDC. This

CE

means that interrupting at the TDC shifts the time-averaged position of the COL forward and makes a nose-up torque through one flapping period,

AC

while interrupting at the BDC shifts the time-averaged position of the COL backward and makes a nose-down torque. Therefore, it can be expected from this result that the pitching angle can be controlled by interrupting the wing motion shortly at the TDC and/or the BDC. Next, we calculate the motion of the wing–body model when the input

34

CR IP T

ACCEPTED MANUSCRIPT

AN US

Figure 15: Results when the wing motion is interrupted shortly at the TDC or the BDC

every flapping period: (a) time variation of the pitching angle of the body and (b) trajectory of the COB. In (b), the initial position of the COB is denoted by (xb , yb ) = (0, 0), and the symbols on the trajectories indicate the position of the COB when the wings are at TDC.

M

intervals Ttop and Tbottom are determined by Eqs. (27) and (28). In this simulation, we set Kpt = 0.03T /1◦ = 1.72 × T , Kit = 0.0018/1◦ = 0.103,

ED

Kdt = T 2 /600/1◦ = 0.0955 × T 2 , Kpb = −Kpt , Kib = −Kit , Kdb = −Kdt and θc = 0◦ as an example. Since Ttop and Tbottom should be positive for θpitch ≤ 0

PT

and θpitch ≥ 0, respectively, from the above preliminary calculations, we use positive gains for Ttop and negative gains for Tbottom . Fig. 16 shows the time

CE

variation of the pitching angle of the body, the trajectory of the COB in the x-y plane, and the time variation of the input intervals Ttop and Tbottom .

AC

Note that ‘Without Control’ in Figs. 16(a) and (b) represents the result for Ttop = Tbottom = 0 and θc = 0◦ . It can be seen from Fig. 16(a) that

with the present control the pitching angle can be controlled in the range of |θpitch | < 10◦ beyond 20 strokes, while without the control it increases

rapidly to 90◦ . Fig. 16(b) shows that the wing–body model with the control 35

ACCEPTED MANUSCRIPT

can go upward continuously. We can see from Fig. 16(c) that the wing motion is interrupted at the BDC much more frequently than at the TDC.

CR IP T

This is because the wing–body model without the control receives a nose-up torque through one flapping period, and consequently the nose-down torque generated by interrupting at the BDC is much more required for maintaining the balance.

Without Control Controlled

60

yb / L W

pitch

[deg]

80

(b)

40 20 0 5

10

15 20 t/T

0.1 0.05 0 0.05 0.1 0.15 0.2 0.25 0.3

30

35

-2 0 2 4 6 8 10 12 14 16 18 xb / LW

Ttop

CE

PT

Ttop / T, Tbottom / T

ED

(c)

25

M

0

5 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 -0.5

AN US

(a) 100

Tbottom 0

5

10

15 20 t/T

25

30

35

AC

Figure 16: Results for the pitching motion control by the COL controller: (a) time variation of the pitching angle of the body, (b) trajectory of the COB and (c) time variation of the input intervals. In (b), initial position of the COB is denoted by (xb , yb ) = (0, 0), and the symbols on the trajectories indicate the position of the COB when the wings are at the TDC. In (c), grey and white backgrounds mark Ttop and Tbottom , respectively.

36

ACCEPTED MANUSCRIPT

Finally, we simulate the response to a sudden large disturbance with the COL controller. A large pitching torque Tzd is added as the disturbance in

CR IP T

the right hand side of Eq. (24) during 10T ≤ t ≤ 10.1T . In this simula-

tion, we set Tzd = 1.5 × ρf Utip L3W for simulating the response to a nose-up

disturbance and Tzd = −1.5 × ρf Utip L3W for simulating the response to a nose-down disturbance. As for tuning parameters for the control, we set

Kpt = −Kpb = 1.72 × T , Kit = Kib = 0, and Kdt = −Kdb = 0.0955 × T 2 . In

AN US

addition, we set θc = 0◦ at first, and we set θc = −10◦ when the vertical velocity of the body at the TDC is smaller than −0.04Utip . Fig. 17 shows

the time variation of the pitching angle of the body and the trajectory of the COB in the x-y plane. It is seen from Fig. 17(a) that the wing–body

80 40

ED

pitch

[deg]

60 20 0 -40 -60

PT

-20

5

CE

0

(b)

yb / LW

(a)

M

model with the COL controller can recover its attitude from the large nose-up

Nose-up disturbance Nose-down disturbance

10 15 20 25 30 35 40 t/T

1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -2 0 2 4 6 8 10 12 14 16 18 20 xb / LW

AC

Figure 17: Results for the pitching motion control by the COL controller with a sudden large disturbance: (a) time variation of the pitching angle of the body and (b) trajectory of the COB. In (b), initial position of the COB is denoted by (xb , yb ) = (0, 0), and the

symbols on the trajectories indicate the position of the COB when the wings are at the TDC.

37

ACCEPTED MANUSCRIPT

and nose-down disturbances during 10T ≤ t ≤ 20T . However, we can see in Fig. 17(b) that the wing–body model largely goes down after the disturbance

CR IP T

due to the disturbed attitude and a long pause of the wing motion during 10T ≤ t ≤ 20T . In this duration, since the vertical velocity of the body at

the TDC is smaller than −0.04Utip , the center of the flapping angle θc is set to −10◦ . Due to the negative center of the flapping angle, the wing–body

model can go upward again as shown in Fig. 17(b). From the results we find

AN US

that the pitching motion with the COL controller works well even for the large disturbance.

5.5. Flight control of actual butterflies and practical flapping air vehicles At the end of this section, we mention the flight control of actual but-

M

terflies and practical flapping air vehicles. Although the present wing–body model is far from actual butterflies and practical flapping air vehicles due to

ED

the simplified wing shape and kinematics and the absence of deformation, the effects of these factors on the flight control might essentially result in the change of the relative position between the COM and the COL. Since

PT

the present COM and COL controllers can appropriately adjust the relative position between the COM and the COL by the PID control scheme, it can

CE

be expected that the concepts of the present controllers should work well and will be useful control methods even for practical flapping air vehicles. Also,

AC

the results for the COL controller suggest a possibility that actual butterflies might adjust the relative position between the COM and the COL by changing the wing kinematics slightly.

38

ACCEPTED MANUSCRIPT

6. Conclusions We have discussed the flight control of the butterfly-like flapping wing–

CR IP T

body model. In this study, we propose two control methods based on the change of the relative position between the center of mass (COM) and the

center of lift (COL), namely the COM controller and the COL controller. In the COM controller, the COM of the wing–body model is moved actively by using a moving weight attached to the body of the wing–body model. In

AN US

the COL controller, the wing motion is interrupted shortly at the top dead

center (TDC) and/or the bottom dead center (BDC) of the flapping motion. The following three kinds of simulations for Re = 300 were performed by using the IB-LBM: (i) the pitching motion control with the COM controller, (ii) the rolling motion control with the COM controller, and (iii) the pitching

M

motion control with the COL controller.

First, we simulated the pitching motion control by actively moving a

ED

weight along the body. The position of the weight is determined by a P control scheme. As a result, we found that the pitching angle is controlled

PT

in the range of |θpitch | < 10◦ , and the wing–body model can go upward against gravity. Second, we simulated the rolling motion control from a

CE

large initial disturbance on the rolling angle by actively moving a weight perpendicular to the body. The position of the weight is determined by a

AC

PI control scheme. As a result, we found that the rolling angle is controlled in the range of |θroll | < 5◦ from the initial disturbance, and the wing–body model can go upward with a small lateral motion. Finally, we simulated the pitching motion control by interrupting the wing motion shortly at the TDC and/or the BDC. The interval of the interruption is determined by a PID 39

ACCEPTED MANUSCRIPT

control scheme. As a result, we found that the pitching angle is controlled in the range of |θpitch | < 10◦ . In addition, it was found that the pitching motion

CR IP T

control with the COL controller works well even for a large disturbance. Controlling the pitching motion and the rolling motion simultaneously is remained in future work. Acknowledgement

AN US

This research used computational resources of the HPCI system provided by ACCMS (Kyoto University) through the HPCI System Research Project (Project ID: hp140025).

Appendix A. Summary of simulations

M

We summarize all cases considered in the present study in Table A.1.

AC

CE

PT

ED

More detailed conditions and results are found in the corresponding sections.

40

CR IP T

ACCEPTED MANUSCRIPT

Table A.1: Summary of simulations. ‘Without’ represents the case without controller,

‘Prelim.’ represents the preliminary calculations, and ‘Disturbed’ represents the case for

Case

Body motion

Wing motion

Control

Result

5.1

Without

Eqs. (22)–(24)

Eqs. (1), (2)

N/A

Fig. 8

Prelim.

Eqs. (8)–(10)

Eqs. (1), (2)

` is fixed

Fig. 9

Control (i)

Eqs. (8)–(10)

Eqs. (1), (2)

Eq. (11)

Fig. 10

Prelim.

Eqs. (16)–(19)

Eqs. (1), (2)

` is fixed

Fig. 11

Control (ii)

Eqs. (16)–(19)

Eqs. (1), (2)

Eq. (20)

Fig. 12

Prelim. 1

Body is fixed

Eqs. (1), (2)

N/A

Fig. 13

Prelim. 2

Eqs. (22)–(24)

Eqs. (25), (26)

θc is fixed

Fig. 14

Eqs. (22)–(24)

Eqs. (25), (26)

Ttop/bottom is fixed

Fig. 15

Control (iii) Eqs. (22)–(24)

Eqs. (25), (26)

Eqs. (27), (28)

Fig. 16

Disturbed

Eqs. (22)–(24)

Eqs. (25), (26)

Eqs. (27), (28)

Fig. 17

5.2

ED

5.3

M

Section

PT

AN US

the control (iii) with a sudden large disturbance.

Prelim. 3

AC

CE

5.4

41

ACCEPTED MANUSCRIPT

References [1] K. Y. Ma, P. Chirarattananon, S. B. Fuller, R. J. Wood, Controlled

CR IP T

flight of a biologically inspired, insect-scale robot, Science 340 (2013) 603–607.

[2] L. Ristroph, S. Childress, Stable hovering of a jellyfish-like flying machine, J. R. Soc. Interface 11 (2014) 20130992 (7pp).

AN US

[3] M. Sun, Insect flight dynamics: Stability and control, Rev. Mod. Phys. 86 (2014) 615–646.

[4] G. K. Taylor, A. L. R. Thomas, Dynamic flight stability in the desert locust Schistocerca gregaria, J. Exp. Biol. 206 (2003) 2803–2829.

M

[5] M. Sun, Y. Xiong, Dynamic flight stability of a hovering bumblebee, J.

ED

Exp. Biol. 208 (2005) 447–459.

[6] M. Sun, J. K. Wang, Y. Xiong, Dynamic flight stability of hovering

PT

insects, Acta Mech. Sin. 23 (2007) 231–246. [7] Y. L. Zhang, M. Sun, Dynamic flight stability of a hovering model insect:

CE

lateral motion, Acta Mech. Sin. 26 (2010) 175–190. [8] I. Faruque, J. S. Humbert, Dipteran insect flight dynamics. Part 1: Lon-

AC

gitudinal motion about hover, J. Theor. Biol. 264 (2010) 538–552.

[9] I. Faruque, J. S. Humbert, Dipteran insect flight dynamics. Part 2: Lateral-directional motion about hover, J. Theor. Biol. 265 (2010) 306– 313.

42

ACCEPTED MANUSCRIPT

[10] M. Sun, J. K. Wang, Flight stabilization control of a hovering model insect, J. Exp. Biol. 210 (2007) 2714–2722.

CR IP T

[11] Y. L. Zhang, M. Sun, Stabilization control of a hovering model insect: lateral motion, Acta Mech. Sin. 27 (2011) 823–832.

[12] B. Cheng, X. Deng, T. L. Hedrick, The mechanics and control of pitching

manoeuvres in a freely flying hawkmoth (Manduca sexta), J. Exp. Biol.

AN US

214 (2011) 4092–4106.

[13] L. Ristroph, A. J. Bergou, G. Ristroph, K. Coumes, G. J. Berman, J. Guckenheimer, Z. J. Wang, I. Cohen, Discovering the flight autostabilizer of fruit flies by inducing aerial stumbles, Proc. Natl. Acad. Sci.

M

U.S.A. 107 (2010) 4820–4824.

[14] L. Ristroph, G. Ristroph, S. Morozova, A. J. Bergou, S. Chang, J. Guck-

ED

enheimer, Z. J. Wang, I. Cohen, Active and passive stabilization of body pitch in insect flight, J. R. Soc. Interface 10 (2013) 20130237 (13pp).

PT

[15] H. Liu, C. P. Ellington, K. Kawachi, C. van den Berg, A. P. Willmott, A computational fluid dynamics study of hawkmoth hovering, J. Exp.

CE

Biol. 201 (1998) 461–477. [16] H. Liu, K. Kawachi, A numerical study of insect flight, J. Comput. Phys.

AC

146 (1998) 124–156.

[17] H. Aono, F. Liang, H. Luo, Near- and far-field aerodynamics in insect hovering flight: an integrated computational study, J. Exp. Biol. 211 (2008) 239–257. 43

ACCEPTED MANUSCRIPT

[18] H. Liu, Integrated modeling of insect flight: From morphology, kinematics to aerodynamics, J. Comput. Phys. 228 (2009) 439–459.

CR IP T

[19] N. Gao, H. Aono, H. Liu, Perturbation analysis of 6DoF flight dynamics and passive dynamic stability of hovering fruit fly Drosophila melanogaster, J. Theor. Biol. 270 (2011) 98–111.

[20] N. Yokoyama, K. Senda, M. Iima, N. Hirai, Aerodynamic forces and

25 (2013) 021902 (24pp).

AN US

vortical structures in flapping butterfly’s forward flight, Physics of Fluids

[21] D. Wu, K. S. Yeo, T. T. Lim, A numerical study on the free hovering flight of a model insect at low Reynolds number, Comput. Fluids 103 (2014) 234–261.

M

[22] K. Ota, K. Suzuki, T. Inamuro, Lift generation by a two-dimensional symmetric flapping wing: immersed boundary–lattice Boltzmann simu-

ED

lations, Fluid Dyn. Res. 44 (2012) 045504 (27pp).

PT

[23] Y. Kimura, K. Suzuki, T. Inamuro, Flight simulations of a twodimensional flapping wing by the IB-LBM, Int. J. Mod. Phys. C 25

CE

(2014) 1340020 (8pp). [24] K. Minami, K. Suzuki, T. Inamuro, Free flight simulations of a

AC

dragonfly-like flapping wing–body model using the immersed boundary– lattice Boltzmann method, Fluid Dyn. Res. 47 (2015) 015505 (17pp).

[25] K. Suzuki, K. Minami, T. Inamuro, Lift and thrust generation by a butterfly-like flapping wing–body model: immersed boundary–lattice Boltzmann simulations, J. Fluid Mech. 767 (2015) 659–695. 44

ACCEPTED MANUSCRIPT

[26] K. Suzuki, T. Inamuro, Effect of internal mass in the simulation of a moving body by the immersed boundary method, Comput. Fluids 49

CR IP T

(2011) 173–187. [27] C. P. Ellington, The novel aerodynamics of insect flight: applications to micro-air vehicles, J. Exp. Biol. 202 (1999) 3439–3448.

[28] L. Ristroph, A. J. Bergou, J. Guckenheimer, Z. J. Wang, I. Cohen,

AN US

Paddling mode of forward flight in insects, Phys. Rev. Lett. 106 (2011) 178103 (4pp).

[29] W. Shyy, Y. Lian, J. Tang, D. Viieru, H. Liu, Aerodynamics of low Reynolds number flyers, Cambridge University Press, New York, 2008.

M

[30] R. Dudley, Biomechanics of flight in neotropical butterflies: morphometrics and kinematics, J. Exp. Biol. 150 (1990) 37–53.

ED

[31] T. Inamuro, Lattice Boltzmann methods for moving boundary flows,

AC

CE

PT

Fluid Dyn. Res. 44 (2012) 024001 (21pp).

45