隨著航空技術的發展日漸迅速，拍撲翼微飛行器已成為目前甚受歡迎的研究主題，一般而言，拍撲翼飛行器具有重量輕及體積小之特點，其小尺度的特性特別適合應用在軍事任務等較隱密用途，此特性對於情搜救難等更具重要性。另一方面，近年來極端氣候變得越來越顯著，因此本篇論文考慮天氣及環境影響，以估算這些因素如何影響拍撲翼飛行器飛行性能，本研究團隊長期專注於天氣環境領域的研究，分析不同氣候條件下之影響程度，為了瞭解拍撲翼飛行器在大雨環境下的影響，吾人推導出其性能方程式，以量化參數來衡量其飛行性能。欲理解拍撲翼運動的物理行為，吾人承襲以往研究成果，首先驗證拍撲翼二維“figure-eight”運動模式，接著驗證三維蝴蝶外形在無風無雨情況下的性能，並考量此蝴蝶身體與翅膀均擺動之狀況，之後對於不同大氣條件下的三維蝴蝶型運動參數進行數值模擬與比較，最後再以克利金最佳化方法來找出傾盆大雨下如何控制頻率變化以維持定速前飛的運動模式。吾人使用基於有限體積法(FVM)的計算流體力學軟體ANSYS Fluent v16.0來模擬蝴蝶周圍的流場現象，並搭配了編譯的UDF(User-Defined-Function)來控制其蝴蝶的拍撲模式，大雨的二相流模式以Discrete Phase Model (DPM)來建構，比較晴朗以及傾盆大雨天氣下的升阻力係數差異,進而利用Matlab計算並比較其性能差異性，再以最佳化Kriging方法算出在傾盆大雨下，必須提升多少的拍撲頻率以將性能提升到晴朗天氣下的原始狀態，最後並比較二者之能量差異值，進而了解能量供應器必須額外付出多少能量，才能將拍撲翼飛行器拉回至晴朗天氣下之飛行模式。相信在本研究中獲得的量化資訊，將在未來惡劣天氣下的拍撲翼飛行器最佳化設計上提供若干設計參考，並了解欲維持其飛行性能而減少的飛行滯空時間。

As the development of aerospace engineering technology progress more rapidly, now one of the most popular topics is Flapping-Wing Aerial Vehicle. Generally, it has light-weight and small size characteristics, and its small scale is particularly suited for the military mission. On the other hand, the recently extreme weather is becoming worse and worse, thus current work considers the climate and environment effects to measure and calculate how much degradation the environment affects the performance of flapping-wing aircraft. Our research team has studied the impact of weather factors for a long time and collected lots experiences in the analysis of different climatic conditions. In order to realize how much environmental impact does affect, we derive the flight vehicle energy equation to measure its flight performance in terms of a quantitative parameter. To understand the flapping motion physical behavior, the 2-D “figure-eight” flapping airfoil motions are first validated, and then the 3-D butterfly performance parameter under different atmospheric conditions are numerically simulated and compared. Both the butterfly body and wing flapping motions are activated and implemented in current research. Our selection of software is ANSYS FLUENT arranging in pairs with UDF (User-Defined-Function) and DPM (Discrete Phase Model) to fulfill this study. After compared the difference between the weather under no rain condition and heavy rain condition, the code which written in Matlab is taken to calculate its performance. From the difference of performance parameter under clear weather and heavy rain condition, it can show how much flapping frequency it has to be increased to maintain the same flight condition under clear weather. Compared the energy loss between clear weather and heavy rain condition to know how much energy the energy supply device should pay under the heavy rain condition. It is believed that the quantitative information gained in this work will be useful in the later optimal design of flapping wing aerial vehicle under some severe weather situations.

Contents
List of Tables	VIII
List of Figures	IX
Nomenclature	XVI
Chapter 1 Introduction	1
1.1 The History of Aviation	1
1.2 Flapping-Wing and Micro Aerial Vehicle (MAV)	4
1.3 Research Motive	7
Chapter 2 Research Background	9
2.1 Flight in Nature	9
2.2 Environmental Impact	10
2.2-1 Heavy Rain	11
2.3 Aerodynamics of Flapping Wing Aircraft	11
2.4 Leading-Edge Vortex and Dynamic Stall	14
2.5 Clap and Fling Mechanism	17
2.6 Flapping Motion	18
2.7 Kriging Method	19
Chapter 3 Numerical Method	21
3.1 Geometry Model Construction	21
3.2 Preprocessing	22
3.3 Governing Equations	25
3.4 Flow Solver and UDF	25
3.5 Dynamic Grid	27
3.6 Discrete Phase Model	28
3.7 Heavy Rain Modeling	29
3.8 Grid Generation	32
3.9 Performance of Flapping Wing	34
3.9-1 3-D F-factor	34
3.9-2 Rainfall Droplet Impact	40
3.10 Kriging Method	42
Chapter 4 Numerical Results	46
4.1 The verification of 2-D simulation	46
4.2 The verification of 3-D simulation	48
4.3 The simulation under heavy rain condition	51
4.4 The energy performance of 3-D butterfly	60
4.4-1 One degree of freedom	61
4.4-2 Two degrees of freedom	67
4.5 Optimization results	74
Chapter 5 Conclusions	97
References	99

List of Tables
Table 1 Parameters of butterfly Morphopeleides mode [14]	24
Table 2 Particle properties in the case of LWC=29g/m3	31
Table 3 The unstructured grids number [12]	34
Table 4 Moment of inertia of wing and body [14]	39
Table 5 Angular velocity of wing	39
Table 6 The parameters of rain droplet	41
Table 7 Classification of the different grids [12] 49
Table 8 Average lift and drag coefficients in different cases	54
Table 9 Cases of different angle of attack	62
Table 10 Average contribution to rate of specific energy of the fifth cycle in case 1	64
Table 11Average contribution to rate of specific energy of the fifth cycle in case 2	67
Table 12 Average contribution to rate of specific energy in different cases	68
Table 13 Drag and lift coefficients compared with different cases of adjusting factors.	75
Table 14 The difference of lift coefficient between Kriging method and CFD calculation	85
Table 15 Average contribution to the rate of specific energy under the condition of frequency=16.52Hz, amplitude=80deg, and LWC= 29g/m3	86
Table 16 The difference of lift coefficient between Kriging method and CFD calculation under LWC=9g/m3	91

List of Figures
Figure 1 A mechanical wing device – ca. 1485 [1]	2
Figure 2 Models in 1870 and 1871[2]	3
Figure 3 University of Toronto flapping-wing vehicle (1991) [3]	3
Figure 4 Wing span vs. body mass, the gray area is the ideal area for MAV [4]	6
Figure 5 Frequency vs. body mass, the gray part is the area for MAV [4]	6
Figure 6 Reynolds number range for flight vehicles [5]	7
Figure 7 “Figure-eight” flapping motion [6]	10
Figure 8 von Karman vortex street [6]	12
Figure 9 Reverse von Karman vortex street [6]	13
Figure 10 In soaring flight, birds use both the updraft thermals and orographic lifting to maintain or gain height and save energy [6]	14
Figure 11 Flapper flow visualization with smoke released from the leading edge of the wing [2]	16
Figure 12 The LEV generation process [9]	17
Figure 13 Clap and fling mechanism [9]	18
Figure 14 The positions of a wing element in one period [10]	19
Figure 15 The size of various parts of the butterfly [12]	21
Figure 16 The 3-D butterfly flapping angle	23
Figure 17 The 3-D butterfly pitching angle	23
Figure 18 The range of 3-D butterfly flapping angle [12]	24
Figure 19 The computational domain [12]  32
Figure 20 Grids and calculating domain of butterfly model [12]33
Figure 21 Grids around the butterfly model in the internal domain [12]	33
Figure 22 Correlation with P and θ [11]	44
Figure 23 Stroke plane angle β and phase angle of translation of wing [6]	47
Figure 24 Lift coefficient of 2-D butterfly in first ten cycles	47
Figure 25 Drag coefficient of 2-D butterfly in first ten cycles	48
Figure 26 Grid convergences in different grid types [12]	50
Figure 27 Lift and drag coefficients vs. time with average CL=0.562, CD=0.0046 in condition of frequency=9Hz, amplitude=80deg, and LWC=0g/m3	51
Figure 28 Particles trace in time of butterfly model from the local view	53
Figure 29 Particles trace in time of butterfly model along the X-axis	53
Figure 30 Lift and drag coefficients vs. time in the condition of frequency=9Hz, amplitude=80deg, and LWC=29g/m3	54
Figure 31 Pressure contour compared in the upper surface with the case of LWC=0g/m3 on the left and the case of LWC=29g/m3 on the right at 0.25T	55
Figure 32 Pressure contour compared in the lower surface with the case of LWC=0g/m3 on the left and the case of LWC=29g/m3 on the right at 0.25T	55
Figure 33 Vorticity contour compared in the upper surface with the case of LWC=0g/m3 on the left and the case of LWC=29g/m3 on the right at 0.25T	56
Figure 34 Vorticity contour compared in the lower surface with the case of LWC=0g/m3 on the left and the case of LWC=29g/m3 on the right at 0.25T	56
Figure 35 Pressure contour compared in the upper surface with the case of LWC=0g/m3 on the left and the case of LWC=29g/m3 on the right at 0.75T	57
Figure 36 Pressure contour compared in the lower surface with the case of LWC=0g/m3 on the left and the case of LWC=29g/m3 on the right at 0.75T	57
Figure 37 Vorticity contour compared in the upper surface with the case of LWC=0g/m3 on the left and the case of LWC=29g/m3 on the right at 0.75T	58
Figure 38 Vorticity contour compared in the lower surface with the case of LWC=0g/m3 on the left and the case of LWC=29g/m3 on the right at 0.75T	58
Figure 39 Vortex core region by lambda-2 criterion on the case of LWC=0g/m3 with pressure contour at T=0	59
Figure 40 Vortex core region by lambda-2 criterion on the case of LWC=29g/m3 with pressure contour at T=0	59
Figure 41 Lift and drag coefficients vs. flapping period in case 1, with average CL=-0.92966 and average CD=-0.00031	62
Figure 42 Rate of climb vs. flapping period in case 1	63
Figure 43 Rate of specific energy vs. flapping period in case 1	63
Figure 44 Lift and drag coefficients vs. flapping period in case 1, with average CL=-0.63286 and average CD=-0.00239	65
Figure 45 Rate of climb vs. flapping period in case 2	66
Figure 46 Rate of specific energy vs. flapping period in case 2	66
Figure 47 Rate of climb vs. time for the case of LWC=0g/m3	69
Figure 48 Rate of specific angular kinetic energy of wing for the case of LWC=0g/m3	69
Figure 49 Rate of specific angular kinetic energy of body for the case of LWC=0g/m3	70
Figure 50 Rate of specific energy of raindrops for the case of	70
Figure 51 Rate of specific energy vs. time for the case of LWC=0g/m3	71
Figure 52 Rate of climb vs. time for the case of LWC=29g/m3	71
Figure 53 Rate of specific angular kinetic energy of wing for the case of LWC=29g/m3	72
Figure 54 Rate of specific angular kinetic energy of body for the case of LWC=29g/m3	72
Figure 55 Rate of specific energy of raindrops for the case of LWC= 29g/m3	73
Figure 56 Rate of specific energy vs. time for the case of LWC=29g/m3	73
Figure 57 Lift and drag coefficients vs. time under the condition of frequency=10Hz, amplitude=95deg, and LWC=0g/m3	76
Figure 58 Lift and drag coefficients vs. time under the condition of frequency=10Hz, amplitude=110deg, and LWC=0g/m3	76
Figure 59 Lift and drag coefficients vs. time under the condition of frequency=10Hz, amplitude=95deg, and LWC=29g/m3	77
Figure 60 Lift and drag coefficients vs. time under the condition of frequency=10Hz, amplitude=80deg, and LWC=29g/m3	77
Figure 61 Lift and drag coefficients vs. time under the condition of frequency=11Hz, amplitude=80deg, and LWC=29g/m3	78
Figure 62 Lift and drag coefficients vs. time under the condition of frequency=12Hz, amplitude=80deg, and LWC=29g/m3	78
Figure 63 Lift and drag coefficients vs. time under the condition of frequency=13Hz, amplitude=80deg, and LWC=29 g/m3	79
Figure 64 Lift and drag coefficients vs. time under the condition of frequency=14Hz, amplitude=80deg, and LWC=29g/m3	79
Figure 65 Lift and drag coefficients vs. time under the condition of frequency=15Hz, amplitude=80deg, and LWC=29g/m3	80
Figure 66 Lift and drag coefficients vs. time under the condition of frequency=16Hz, amplitude=80deg, and LWC=29g/m3	80
Figure 67 Lift and drag coefficients vs. time under the condition of frequency=17Hz, amplitude=80deg, and LWC=29g/m3	81
Figure 68 Lift and drag coefficients vs. time under the condition of frequency=18Hz, amplitude=80deg, and LWC=29g/m3	81
Figure 69 Lift and drag coefficients vs. time under the condition of frequency=19Hz, amplitude=80deg, and LWC=29g/m3	82
Figure 70 Lift and drag coefficients vs. time under the condition of frequency=20Hz, amplitude=80deg, and LWC=29g/m3	82
Figure 71 Lift coefficient vs. varying frequencies by Kriging method	85
Figure 72 Lift and drag coefficient vs. time under the condition of frequency=16.52Hz, amplitude=80deg, and LWC=29g/m3	86
Figure 73 Rate of climb vs. time under the condition of frequency=16.52Hz, amplitude=80deg, and LWC=29g/m3	87
Figure 74 Rate of specific angular kinetic energy of wing under the condition of frequency=16.52Hz, amplitude=80deg, and LWC=29g/m3	87
Figure 75 Rate of specific angular kinetic energy of body under the condition of frequency=16.52Hz, amplitude=80degree, and LWC= 29g/m3	88
Figure 76 Rate of specific energy of raindrops under the condition of frequency=16.52Hz, amplitude=80degree, and LWC=29g/m3	88
Figure 77 Rate of specific energy under the condition of frequency=16.52Hz, amplitude=80degree, and LWC=29g/m3	89
Figure 78 Lift coefficient vs. varying frequencies by Kriging method under LWC=9g/m3	91
Figure 79 Lift and drag coefficient vs. time under the condition of frequency=10.83Hz, amplitude=80deg, and LWC=9g/m3	92
Figure 80 Rate of climb vs. time under the condition of frequency=10.83Hz, amplitude=80deg, and LWC=9g/m3	92
Figure 81 Rate of specific angular kinetic energy of wing under the condition of frequency=10.83Hz, amplitude=80deg, and LWC=9g/m3	93
Figure 82 Rate of specific angular kinetic energy of body under the condition of frequency=10.83Hz, amplitude=80deg, and LWC=9g/m3	93
Figure 83 Rate of specific energy of raindrops under the condition of frequency=10.83Hz, amplitude=80deg, and LWC=9g/m3	94
Figure 84 Rate of specific energy under the condition of frequency=10.83Hz, amplitude=80deg, and LWC=9g/m3	94
Figure 85 Specific energy vs. time under different weather and frequency conditions	96

[1]	Leonardo Da Vinci-His Flying Machine, Online Available:
http://www.angelfire.com/eletronic/awakening101/leonardo.html.
[2]	Ellington, C.P., van den Berg, C., Willmott, A.P., and Thomas, A.L. R.,“ Leading Edge Vortices in Insect Flight,” Nature, Vol. 384, 1996, pp. 626-630.
[3]	Kang, Y.H., Ma, Z.H. and Lee, W.T., “Development and Drive Mechanism of Ornithopter,” Journal of Engineering and Technology, Vol. 8, No. 4, 2011, pp. 623-641.
[4]	Shyy, W., Berg M. and Ljungqvist D., “Flapping and Flexible Wings for Biological and Micro Air Vehicles,” Progress in Aerospace Sciences, Vol. 36, 1999, pp.455-505.
[5]	Muller, J.T., and DeLaurier, J.D., “An Overview of Micro Air Vehicle Aerodynamics,” Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications, Vol. 195, 2001, pp. 1-9.
[6]	Aerodynamics of Flapping Flight, Online Available: http://www.dicat.unige.it/guerrero/thesischaps/chapter2.pdf
[7]	Rhode, R.V., “Some Effects of Rainfall on Flight of Airplanes and on Instrument Indications,” NACA TN 903, 1941.
[8]	Theodore von Karman, Aerodynamics: Selected Topics in the Light of Their Historical Development, Dover Publications, 2004.
[9]	Sane, S.P., “The Aerodynamics of Insect Flight,” The Journal of Experimental Biology, Vol. 206, August 2003, pp. 4191-4208.
[10]	Wang, J.Z., “Computation of Insect Hovering,” Mathematical Methods in the Applied Sciences, Vol. 24, 2001, pp. 1515-1521.
[11]	Alexander, I.J. Forrester, S. András, and J. K. Andy, Engineering Design via Surrogate Modeling: A Practical Guide, WILEY, 2008.
[12]	Cheng, Y.T., “On the Simulation of 3-D Flapping-Wing Insect’s Flight Performance under Abnormal Atmospheric Conditions,” MS thesis, Tamkang University, 2014.
[13]	Dudley, R., “Biomechanics of Flight in Neotropical Butterflies: Morphometrics Kinematics,” The Journal of Experimental Biology, Vol. 150, 1990, pp. 37-53.
[14]	Huang, H. and Sun, M, “Forward Flight of a Model Butterfly: Simulation by Equations of Motion Coupled with the Navier–Stokes Equations,” Acta Mechanica Sinica, Vol. 28, 2012, pp. 1590-1601.
[15]	Dunham, R.E. Jr., “The Potential Influence of Rain on Airfoil Performance,” Von Karman Institute for Fluid Dynamics, 1987.
[16]	Markowitz, A.M., “Raindrop Size Distribution Expression,” Journal of Applied Meteorology, Vol. 15, 1976, pp. 1029-1031.
[17]	NASA, “The Hazard Index: Langley’s ‘F-factor’https://www.nasa.gov/centers/langley/news/factsheets/F-Factor.html”
[18]	Wang, R.L., “Numerical Simulation of Flapping Airfoil in Gusty Environment,” MS thesis, Tamkang University, 2012.
[19]	Wan, T., and Y.S. Chen, “On the Optimization of Blended Wing Body Aircraft Configuration via the Surrogate Modeling Method,” 29th International Congress of the Aeronautical Sciences, St. Petersburg, Russia, Sept. 7-12, 2014.