進行蝴蝶模擬前進飛行時非均勻流對於其氣動力改變分析，阻力係數部分，正攻角下拍行程係數值為正、上拍行程為負，負攻角下拍行程係數值為負、上拍行程為正；升力係數部分，下拍行程攻角為正時，有助於提高升力，上拍行程攻角為負時，升力係數有最小值；俯仰力矩係數部分，在下拍行程時會造成一個使蝴蝶抬頭的力，當在下拍行程一個正的攻角，會使蝴蝶擁有更大抬頭的力，於上拍行程時會造成蝴蝶一個抬尾的力，當遇到負攻角時會使抬尾的力增大。
　　身體主動俯仰動作的模擬部分，升力與俯仰力矩與身體兩個俯仰角度之間差別不大，只對阻力係數變化有關，對此推測身體俯仰對於向前飛行時的推力有幫助。而流體經過身體主動俯仰飛行中的蝴蝶時，身體俯仰會造成翅膀上的壓力分佈有所改變。

Aerodynamic simulation of Butterfly’s flapping wings in non-uniform flows is performed. In the prediction of the aerodynamic forces, the drag coefficient is shown to be positive when the angle of attack is positive during the downstroke. However, the drag coefficient is shown to be negative during the upstroke. Also, the drag coefficient is negative when the angle of attack is negative during the downstroke., but the drag coefficient is positive during the upstroke ; For lift coefficient, the lift is enhanced during the up-stroke and reaches its minimum when angle of attack is positive. For moment coefficient, the pitching up force in the downstroke is produced. If angle of attack is positive, the pitching moment is enforced. However, the moment of the pitching-down appears and is enhanced when the angle of attack is reduced in the upstroke.
　　In our simulation Butterfly pitching body in forward flight, the distributions of the lift and moment coefficients are still same in the case without the body motion of pitching except the distribution of the drag coefficient. This combined the body motion of pitching with flapping motion is shown to strengthen the thrust during the forward flight of the butterfly.

目錄
目錄	i
表目錄	ii
圖目錄	iii
第一章 前言	1
1.1研究背景	1
1.2文獻回顧	5
1.3研究目的	8
1.4論文架構	9
第二章 物理模型與數值方法	10
2.1  統御方程式	10
2.2  邊界條件	10
2.3  蝴蝶拍撲模型幾何尺寸	11
2.4  運動模式及參數定義	13
2.5  計算區域與網格設置	15
2.6  數值方法與計算流程	16
第三章 模擬結果與討論	18
3.1  燕尾蝶翼於非均勻流場下飛行模擬	18
3.1.1 非均勻流頻率fnon=1Hz	18
3.1.2 非均勻流頻率fnon=2Hz	19
3.1.3 非均勻流頻率fnon=3Hz	20
3.1.4 非均勻流頻率fnon=3.5Hz	20
3.1.5 非均勻流頻率fnon=4Hz	21
3.1.6 非均勻流頻率fnon=5Hz	21
3.1.7 非均勻流頻率fnon=6Hz	22
3.1.8 非均勻流頻率fnon=7Hz	23
3.1.9 比較各非均勻流不同頻率間壓力梯圖及渦流圖關係	24
3.2  燕尾蝶翼身體俯仰動作之飛行模擬	25
3.3 均勻流與非均勻流和身體俯仰之關係比較	26
第四章 結論與未來展望	28
參考文獻	30


 
表目錄
表 4-1 非均勻流升阻力及俯仰力矩於單一拍撲過程比較圖。	32
表 4-2非均勻流渦度於單一拍撲過程比較圖。	32
 
圖目錄
圖 1-1 各科蝴蝶比較圖(a)弄蝶科: Euschemon rafflesia (b)鳳蝶科: 鳳蝶科Teinopalpus imperialis (c)粉蝶科: Pareronia anais (d)灰蝶科: Lycaeides argyrognomon (e)蛺蝶科: Diaethria eluina。	33
圖2 2 本研究所探討的蝴蝶品種Graphium policenes。	34
圖2 3 H.Park et al. [7]文章中的Graphium policenes蝴蝶外型。	34
圖2 4 依照H.Park ea al.[7]文章裡的蝴蝶輪廓略以修改成的三維蝴蝶模型及各部位尺寸示意圖(上圖為俯視圖、下圖為側視圖)。	35
圖2 5 本研究蝴蝶模型身體各部位尺寸示意圖	35
圖2 6 蝴蝶模型翅膀與胸部連接處示意圖(左圖為俯視圖、又圖為前視圖)，左圖中點線處為翅膀與身體胸部連接處，θ=35゜為下拍行程的起始位置，θ=-35゜為上拍行程的起始位置。	36
圖2 7 本研究中蝴蝶模型翅膀拍撲角度ϕ與攻角α關係示意圖。	36
圖2 8 C.R.Betts et al. [1] 研究中所採用的Troides rhadamantus 蝴蝶翅膀輪廓外型。	37
圖2 9 蝴蝶質量與拍撲頻率關係圖	37
圖2 10 模擬計算區域與蝴蝶模型位置關係圖	38
圖2 11 模擬計算區域與網格設置示意圖	38
圖3-12 非均勻流在不同頻率週期下與蝴蝶拍撲週期對應之攻角比較圖，(a)流場頻率1Hz、(b)流場頻率2Hz、(c)流場頻率3Hz、(d)流場頻率3.5Hz、(e)流場頻率4Hz、(f)流場頻率5Hz、(g)流場頻率6Hz、(h)流場頻率7Hz。	40
圖3 13 非均勻流頻率為1Hz，攻角10°~20°於7次拍撲週期內CD、CL、CM之比較圖，(i)為CD、(ii)為CL、(iii)為CM 。	41
圖3 14非均勻流頻率為2Hz，攻角10°~20°於7次拍撲週期內CD、CL、CM之比較圖，(i)為CD、(ii)為CL、(iii)為CM 。	42
圖3 15非均勻流頻率為3Hz，攻角10°~20°於7次拍撲週期內CD、CL、CM之比較圖，(i)為CD、(ii)為CL、(iii)為CM 。	43
圖3 16非均勻流頻率為3.5Hz，攻角10°~20°於7次拍撲週期內CD、CL、CM之比較圖，(i)為CD、(ii)為CL、(iii)為CM 。	44
圖3 17非均勻流頻率為4Hz，攻角10°~20°於7次拍撲週期內CD、CL、CM之比較圖，(i)為CD、(ii)為CL、(iii)為CM 。	45
圖3 18非均勻流頻率為5Hz，攻角10°~20°於7次拍撲週期內CD、CL、CM之比較圖，(i)為CD、(ii)為CL、(iii)為CM 。	46
圖3 19非均勻流頻率為6Hz，攻角10°~20°於7次拍撲週期內CD、CL、CM之比較圖，(i)為CD、(ii)為CL、(iii)為CM 。	47
圖3 20非均勻流頻率為7Hz，攻角10°~20°於7次拍撲週期內CD、CL、CM之比較圖，(i)為CD、(ii)為CL、(iii)為CM 。	48
圖3-21非均勻流頻率1Hz、最大攻角為20°，上翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	50
圖3-22 非均勻流頻率1Hz、最大攻角為20°，下翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	52
圖3-23 非均勻流頻率1Hz、最大攻角為20°，X方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	54
圖3-24 非均勻流頻率1Hz、最大攻角為20°，Y方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	56
圖3-25 非均勻流頻率1Hz、最大攻角為20°，Z方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	58
圖3-26非均勻流頻率2Hz、最大攻角為20°，上翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	60
圖3-27非均勻流頻率2Hz、最大攻角為20°，下翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	62
圖3-28非均勻流頻率2Hz、最大攻角為20°，X方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	64
圖3-29非均勻流頻率2Hz、最大攻角為20°，Y方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	66
圖3-30非均勻流頻率2Hz、最大攻角為20°，Z方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	68
圖3-31非均勻流頻率3Hz、最大攻角為20°，上翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	70
圖3-32非均勻流頻率3Hz、最大攻角為20°，下翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	72
圖3-33非均勻流頻率3Hz、最大攻角為20°，X方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	74
圖3-34非均勻流頻率3Hz、最大攻角為20°，Y方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	76
圖3-35非均勻流頻率3Hz、最大攻角為20°，Z方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	78
圖3-36 非均勻流頻率3.5Hz、最大攻角為20°，上翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	80
圖3-37非均勻流頻率3.5Hz、最大攻角為20°，下翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	82
圖3-38 非均勻流頻率3.5Hz、最大攻角為20°，X方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	84
圖3-39 非均勻流頻率3.5Hz、最大攻角為20°，Y方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	86
圖3-40非均勻流頻率3.5Hz、最大攻角為20°，Z方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	88
圖3-41 非均勻流頻率4Hz、最大攻角為20°，上翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	90
圖3-42非均勻流頻率4Hz、最大攻角為20°，下翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	92
圖3-43非均勻流頻率4Hz、最大攻角為20°，X方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	94
圖3-44非均勻流頻率4Hz、最大攻角為20°，Y方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	96
圖3-45非均勻流頻率4Hz、最大攻角為20°，Z方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	98
圖3-46非均勻流頻率5Hz、最大攻角為20°，上翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	100
圖3-47非均勻流頻率5Hz、最大攻角為20°，下翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	102
圖3-48非均勻流頻率5Hz、最大攻角為20°，X方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	104
圖3-49非均勻流頻率5Hz、最大攻角為20°，Y方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	106
圖3-50非均勻流頻率5Hz、最大攻角為20°，Z方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	108
圖3-51非均勻流頻率6Hz、最大攻角為20°，上翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	110
圖3-52非均勻流頻率6Hz、最大攻角為20°，下翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	112
圖3-53非均勻流頻率6Hz、最大攻角為20°，X方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	114
圖3-54非均勻流頻率6Hz、最大攻角為20°，Y方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	116
圖3-55非均勻流頻率6Hz、最大攻角為20°，Z方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	118
圖3-56非均勻流頻率7Hz、最大攻角為20°，上翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	120
圖3-57非均勻流頻率7Hz、最大攻角為20°，下翼面壓力梯圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	122
圖3-58非均勻流頻率7Hz、最大攻角為20°，X方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	124
圖3-59非均勻流頻率7Hz、最大攻角為20°，Y方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	126
圖3-60非均勻流頻率7Hz、最大攻角為20°，Z方向渦度圖，(a)第一拍撲週期、(b)第二拍撲週期、(c)第三拍撲週期、(d)第四拍撲週期、(e)第五拍撲週期、(f)第六拍撲週期、(g)第七拍撲週期。	128
圖 3-61 俯仰動作於攻角10度、20度、30度下的第三拍撲週期CD、CL、CM之比較圖，(i)為CD、(ii)為CL、(iii)為CM 。	129
圖 3-62 俯仰動作於第三拍撲週期，最大攻角為20°，上翼面壓力梯圖。	130
圖 3-63 俯仰動作於第三拍撲週期，最大攻角為20°，下翼面壓力梯圖。	130
圖 3-64 俯仰動作於第三拍撲週期，最大攻角為20°，X方向渦度圖。	130
圖 3-65 俯仰動作於第三拍撲週期，最大攻角為20°，Y方向渦度圖。	131
圖 3-66 俯仰動作於第三拍撲週期，最大攻角為20°，Z方向渦度圖。	131
圖 3-67 均勻流與非均勻流和身體俯仰於攻角10度下第三拍撲週期CD、CL、CM之比較圖，(i)為CD、(ii)為CL、(iii)為CM 。	132

參考文獻
[1]   C.R. Betts , R.J.Wootom.Wing Shape and Behavior In Butterflies: A Preliminary Analysis. J.Exp.Biol.138, 271-288, 1988.
[2]   A.K.Brodsky .Vortex Formation in the Tethered Flight of the Peacock Butterfly Inachis io L. (Lepidoptera, Nymphalidae) and some Aspects of Insect Flight Evolution. J Exp Biol 161, 77-95,1991
[3]   S. Sunada. Preformance of A Butterfly In Take-Off Flight.J.exp.Biol.183,249-277,1993.
[4]   C.P. Ellington, C. Berg, A. P. Willmott and A. R. Thomas. Leading-Edge Vorticesin Insect Flight. Nature 384, 626 – 630,1996
[5]   K. Kawachi and H. Liu. A Numerical Study of Insect Flight. Journal of Computational Physics 146, 124–156 (1998)
[6]   K.Senda , T.Tsuyoshi and M.Sawamoto, Measurement and Numerical Simulation of a Flapping Butterfly. Proceeding of 2nd International Symposium on Adaptive Motion of Animals and Machines,Kyoto,March.4-8,2003
[7]   H. Park ,K. Bae , B. Lee ,W.P. Jeon and H. Choi. Aerodynamic Performance of a Gliding Swallowtail Butterfly Wing Model. Exp Mechanics 50,1313-1321,2010.
[8]   H. Tanaka and I. Shimoyama. Forward Flight of Swallowtail Butterfly with SimpleFlapping Motion Bioinspire. Biomim. 5 026003, 2010
[9]   Hua Huang and Mao Sun. Forward flight of a model butterfly: Simulation by equations of motion coupled with the Navier–Stokes equations. Acta Mechanica Sinica 28(6):1590–1601,2012
[10]  張雅涵,碟型拍撲翼之氣動力研究. 國立台灣大學工學院機械工程學系,2012
[11]  N.Yokoyama, K.Senda, M.Iima, and N.Hirai. Aerodynamic forces and vortical structures in flapping butterfly's forward flight. Phys. of Fluids., 25, 021902,2013.
[12]  M.Fuchiwaki, T.Kuroki, K.Tanaka and T.Tababa. Dynamic behavior of the vortex ring formed on a butterfly wing. Exp Fluids,54:1450,2013..
[13]  R. Dudley. Biomechanics of Flight In Neotropical Butterflies : Morphometrics　Kinematics. J. exp. Biol.150 ,37-53, 1990.
[14]  R. I. Issa. Solution of Implicitly Discretized Fluid Flow Equations by Operator　Splitting. J. Comput. Phys., 62:40–65, 1986.