螺旋槳渦格法的初步實(shí)現(xiàn)

螺旋槳渦格法的初步實(shí)現(xiàn)螺旋槳升力面理論有渦分布法，偶極子分布法和加速度勢(shì)法等。經(jīng)過(guò)多年的應(yīng)用實(shí)踐，現(xiàn)在普遍采用的是渦分布法。雖然面元法的發(fā)展比升力面方法在理論上更為完善，尤其是壓力分布的預(yù)報(bào)計(jì)算顯著優(yōu)于后者，但就螺旋槳合力的計(jì)算而言，兩者基本具有同樣的計(jì)算精度，因此升力面理論方法目前仍是螺旋槳理論設(shè)計(jì)所采用的主要工具。隨著計(jì)算機(jī)速度和內(nèi)存的不斷提高和擴(kuò)大，渦格法(VortexLatticeMethod，即VLM)得到越來(lái)越廣泛的應(yīng)用。至今，各學(xué)者所發(fā)表的各種渦升力面理論的論文，就原理而言都是一致的，而主要差別僅在于渦模型包括尾渦處理和數(shù)值方法兩方面。作為螺旋槳理論的初學(xué)者，為了深入理解升力面理論以及應(yīng)用方式，我們將通過(guò)自編程序在計(jì)算機(jī)上進(jìn)行簡(jiǎn)單的渦格法實(shí)現(xiàn)。由于螺旋槳的幾何形狀較為復(fù)雜，而我們的目的在于對(duì)升力面渦格法有初步的認(rèn)識(shí)，因此我們對(duì)螺旋槳作一定程度上的簡(jiǎn)化，既能反映渦格法的計(jì)算特征，又能使計(jì)算相對(duì)簡(jiǎn)便。在對(duì)渦格法理論模型進(jìn)行簡(jiǎn)單的介紹之后，本文將給出渦格法的實(shí)現(xiàn)程序，并對(duì)不同攻角以及不同拱弧面下的計(jì)算結(jié)果進(jìn)行分析比較。一、渦格法理論模型1、螺旋槳幾何特征本文使用一空間曲面(拱弧面，即升力面)代替螺旋槳槳葉，槳葉厚度可由離散的源匯分布來(lái)模擬，但在本文中不予考慮。該拱弧面在xoy平面的投影為展長(zhǎng)2.0,弦長(zhǎng)1.0的矩形平面，在xoz平面上的投影為一剖面拱弧。為了比較不同的拱弧線對(duì)槳葉性能的影響，本文使用三種不同的剖面拱弧線，它們分別為NACAa二0.8剖面拱弧、圓弧以及拋物線拱弧。為了比較不同的最大拱高值對(duì)槳葉性能的影響，本文取兩種不同的最大拱高，分別為fmax=0.02*c和fmax=0.04*c，其中c為弦長(zhǎng)1.0。2、控制方程渦格法的思想是將升力面設(shè)置在拱弧面上，然后在升力面上進(jìn)行渦格劃分，并在每個(gè)渦格上布置馬蹄渦和控制點(diǎn)，而每個(gè)馬蹄渦的強(qiáng)度是未知的，可以通過(guò)在控制點(diǎn)處滿(mǎn)足拱弧面上法向速度為零的邊界條件進(jìn)行求解。本文中渦格為均勻網(wǎng)格劃分，展向劃分m個(gè)渦格，弦向劃分n個(gè)渦格，在每個(gè)渦格中將馬蹄渦的展向渦段布置在每格的1/4弦向格長(zhǎng)位置處，而控制點(diǎn)置于3/4弦向格長(zhǎng)及1/2展向格長(zhǎng)位置，Weissinger證明了這樣做已隱含著庫(kù)塔條件的滿(mǎn)足。對(duì)于弦向渦線則沿渦格線分成直線渦段，亦即說(shuō)整個(gè)弦向渦的幾何形狀用折線代替，尾渦區(qū)渦線另行處理。這樣，對(duì)于每一個(gè)控制點(diǎn)，滿(mǎn)足物面不可穿透條件可得：(v+v)?n=0IJ0IJ其中v和n分別表示第(I,J)個(gè)渦格控制點(diǎn)的誘導(dǎo)速度和法向矢量，v表示流場(chǎng)的入流速I(mǎi)Jij - - - 0度。v為所有的升力面馬蹄渦在第(I,J)個(gè)渦格控制點(diǎn)的誘導(dǎo)速度的疊加，即：?口f -vIJ+K(vIJ+K(uci+1,l—Uc)+(uwi,l i+1,j-uw)]*ri,j i,j其中us,Y(uc-uc),uw-uw分別表示單位渦強(qiáng)的第(i,j)個(gè)馬蹄渦的展向渦段、i,j i+1,1 i,l i+1,ji,jl=j弦向渦段以及尾渦段對(duì)第(I,J個(gè)渦格控制點(diǎn)的誘導(dǎo)速度，可根據(jù)Biot-Savart定理求得，而廠表示第(i,j)個(gè)馬蹄渦的渦強(qiáng)，為待求的未知量。i,j對(duì)所有的控制點(diǎn)滿(mǎn)足方程聯(lián)立求解，即可得每個(gè)渦格內(nèi)馬蹄渦的渦強(qiáng)。3、 葉面區(qū)的渦系模型在葉面區(qū)內(nèi)，展向均勻劃分15個(gè)格，弦向均勻劃分20格，而每一個(gè)馬蹄渦按渦格劃分為展向渦段和弦向渦段。由于NACAa二0.8的剖面拱弧型值為一系列散點(diǎn)坐標(biāo)，所以在求解渦段端點(diǎn)及控制點(diǎn)坐標(biāo)時(shí)需要進(jìn)行插值處理。考慮到最大拱高值相對(duì)于弦長(zhǎng)較小，使用線性插值和三次樣條插值的結(jié)果差別不大，因此選用較為簡(jiǎn)單的線性插值即可。對(duì)于控制點(diǎn)的法向矢量，通過(guò)求取與控制點(diǎn)相鄰兩點(diǎn)所組成的直線斜率，然后再作簡(jiǎn)單的三角處理即可求得。對(duì)于圓弧以及拋物線拱弧，可直接求得各點(diǎn)的坐標(biāo)值和法向矢量。但為了保持程序的統(tǒng)一性，本文依然先求得拱弧線的一系列散點(diǎn)坐標(biāo)，然后按與NACAa=0.8拱弧相類(lèi)似的處理方法求取各點(diǎn)的坐標(biāo)值和法向矢量。4、 尾流區(qū)的渦系模型尾流區(qū)內(nèi)尾渦的處理較為復(fù)雜，且對(duì)計(jì)算結(jié)果有較大的影響。尾流區(qū)是自由空間，故渦的走向須沿著流場(chǎng)內(nèi)該點(diǎn)的速度方向。而流場(chǎng)速度是需問(wèn)題解出后才知道，故尾流區(qū)是事先未知的，這就構(gòu)成了解升力面問(wèn)題的非線性。為了使求解問(wèn)題簡(jiǎn)化，本文事先假定了尾流區(qū)的范圍和形狀。本文將尾流區(qū)分為過(guò)渡尾流區(qū)及遠(yuǎn)尾流區(qū)，以便更好的符合實(shí)際情況。在過(guò)渡尾流區(qū)，尾渦在隨邊處沿剖面拱弧的切線方向泄出，由于受來(lái)流速度的影響，尾渦在過(guò)渡尾流區(qū)結(jié)束的地方與來(lái)流方向一致。過(guò)渡尾流區(qū)的長(zhǎng)度為一倍弦長(zhǎng)。在過(guò)渡尾流區(qū)內(nèi)，尾渦的走向成二次拋物線型。與葉面區(qū)類(lèi)似，過(guò)渡尾流區(qū)內(nèi)的尾渦按弦向共劃分n個(gè)渦段。同時(shí)，考慮到展向兩端的尾渦卷曲，尾渦在隨邊處的下泄方向不僅與剖面拱弧的切線方向有關(guān)，而且與來(lái)流速度相關(guān)，因此本文中尾渦的下泄方向?yàn)椋簐=av+(1—ci)vTOC\o"1-5"\h\zw 0 t其中V，v，v分別表示尾渦下泄方向，來(lái)流速度方向和隨邊處剖面拱弧的切線方向，cw0t 一 f 一與展向坐標(biāo)有關(guān)，其表達(dá)式為：f f fa=2*(y/s—0.5)2其中y為尾渦下泄點(diǎn)的展向坐標(biāo)，s為展長(zhǎng)。在遠(yuǎn)尾流區(qū)，尾渦方向與來(lái)流速度方向一致，為半無(wú)限長(zhǎng)渦線，處理較為簡(jiǎn)單。二、數(shù)值結(jié)果的比較分析

本文計(jì)算了NACAa二0.8剖面拱弧、圓弧以及拋物線等三種不同的拱弧面，不同最大拱高值分別為fmax=0.02*c和fmax=0.04*c（其中c為弦長(zhǎng)1.0）,不同攻角分別為0.0，2.5,5.0和10.0下的渦強(qiáng)分布。為了便于比較，本文中所給結(jié)果僅為槳葉在展向中部的渦強(qiáng)分布，此處的渦強(qiáng)值最大。下面給出部分?jǐn)?shù)值結(jié)果的比較分析，其中離散點(diǎn)標(biāo)識(shí)為數(shù)值計(jì)算結(jié)果，而線條值為樣條擬合結(jié)果，由于渦強(qiáng)在導(dǎo)邊附近變化較為激烈，故擬合值可能存在較大誤差，僅作為參考。1、不同攻角下渦強(qiáng)沿弦向分布的比較圖1所示為當(dāng)剖面拱弧和最大拱高一定時(shí)，不同攻角下渦強(qiáng)沿弦向分布的比較。從圖1中可以看出，當(dāng)攻角從0.0增加到10.0度時(shí)，渦強(qiáng)均增大，且在導(dǎo)邊處變化較為明顯，隨邊處變化較小。當(dāng)fmax=0.02*c時(shí)，攻角為0.0導(dǎo)致導(dǎo)邊附近環(huán)量為負(fù)值；當(dāng)fmax=0.04*c時(shí)，攻角為0.0和2.5均導(dǎo)致導(dǎo)邊附近環(huán)量為負(fù)值，說(shuō)明具有較大拱高值的拱弧線在大攻角入流條件下才能取得較好的升力性能。同時(shí)，NACAa二0.8剖面拱弧在0?0.8倍弦長(zhǎng)范圍內(nèi)渦強(qiáng)變化較為平緩，在0.8~1.0倍弦長(zhǎng)內(nèi)渦強(qiáng)呈直線下降到0，且?guī)焖l件滿(mǎn)足較好；而圓弧和拋物線剖面在整個(gè)弦長(zhǎng)范圍內(nèi)渦強(qiáng)變化較為平緩，但庫(kù)塔條件滿(mǎn)足不是很好。(a)(a)NACAa=0.8剖面，fmax=0.02*c(b)NACAa二0.8剖面，fmax=0.04*ca-100yisioo.*-?■00??00下庚攻甬下潘區(qū)沿蔻向分布的比較■:arc:fmax=0.02-ca-100yisioo.*-?■00??00下庚攻甬下潘區(qū)沿蔻向分布的比較■:arc:fmax=0.02-c）下潘強(qiáng)沿處向分布？H匕較arc:fmax=OO4-c>-002(c)圓弧剖面,(c)圓弧剖面,fmax=0.02*c(d)圓弧剖面，fmax=0.04*c0.1"100iM-a?00不同攻角下?tīng)?wèi)強(qiáng)沿菠向分心的比較parabolafmax=002*c0.1"100iM-a?00不同攻角下?tīng)?wèi)強(qiáng)沿菠向分心的比較parabolafmax=002*c不同攻角下?tīng)?wèi)強(qiáng)沿菠向分專(zhuān)的比較parabolafmax=004*c-002(e)拋物線剖面，fmax=0.02*c(f)拋物線剖面，fmax=0.04*c圖1、不同攻角下渦強(qiáng)沿弦向分布的比較2、不同拱弧下渦強(qiáng)沿弦向分布的比較圖2所示為當(dāng)最大拱高值和攻角一定時(shí)，不同拱弧下渦強(qiáng)沿弦向分布的比較。從圖中可以看到，圓弧剖面和拋物線剖面所得到的渦強(qiáng)分布極其相似，幾乎無(wú)法分辨出兩者的差別可能需要增加插值點(diǎn)數(shù)才能提高精度，顯示出兩者的差別。同時(shí)，NACAa0.8剖面拱弧在0~0.8倍弦長(zhǎng)范圍內(nèi)渦強(qiáng)變化較為平緩，在0.8~1.0倍弦長(zhǎng)內(nèi)渦強(qiáng)呈直線下降到0，且?guī)焖l件滿(mǎn)足較好；而圓弧和拋物線剖面在整個(gè)弦長(zhǎng)范圍內(nèi)渦強(qiáng)都有變化，且?guī)焖l件滿(mǎn)足不是很好。當(dāng)fmax=0.04*c，攻角為2.5度時(shí)，螺旋槳未達(dá)到正常的工作狀態(tài)。此外，圓弧剖面和拋物線剖面所得到的渦強(qiáng)分布擬合曲線出現(xiàn)鋸齒狀，可能是由于求剖面型值時(shí)插值節(jié)點(diǎn)較少導(dǎo)致精度不夠，增加插值節(jié)點(diǎn)可能會(huì)有所改善。(a)fmax=0.02*c,攻角為2.5度 (b)fmax=0.02*c,攻角為5.0度(c)(c)fmax=0.04*c，攻角為2.5度(d)fmax=0.04*c，攻角為5.0度圖2、不同拱弧下渦強(qiáng)沿弦向分布的比較3、不同拱高下渦強(qiáng)沿弦向分布的比較3、不同拱高下渦強(qiáng)沿弦向分布的比較圖3所示為NACAa二0.8剖面在攻角一定時(shí)，不同最大拱高值下渦強(qiáng)沿弦向分布的比較。從圖中可以看到，當(dāng)fmax=0.02*c，攻角為2.5度，或fmax=0.04*c，攻角為5.0度時(shí)，NACAa二0.8剖面得到的渦強(qiáng)分布較為理想。(a)(a)NACAa二0.8剖面，攻角為0.0度(b)NACAa二0.8剖面，攻角為2.5度0.04*fmax?O02"cfmax?OO4*c02ac002060406不同拱&下潘暹檔菠向分布的比較NACA:8=10.0.不同拱玄下?tīng)?wèi)強(qiáng)縉菠向分畝的比較NACA;a=5.0003500250.04*fmax?O02"cfmax?OO4*c02ac002060406不同拱&下潘暹檔菠向分布的比較NACA:8=10.0.不同拱玄下?tīng)?wèi)強(qiáng)縉菠向分畝的比較NACA;a=5.0003500250.020005-0015,fmax=0M*c"(a)(a)NACAa二0.8剖面，攻角為5.0度(b)NACAa二0.8剖面，攻角為10.0度圖3、不同拱高下渦強(qiáng)沿弦向分布的比較三、渦格法程序的實(shí)現(xiàn)及數(shù)值結(jié)果的后處理本文中渦格法程序使用C語(yǔ)言編寫(xiě)，起初的編寫(xiě)環(huán)境為linux操作系統(tǒng)，編譯工具為gcc。后來(lái)為了使用Matlab作圖，于是改用Windows操作系統(tǒng)的MicrosoftVisualC++編譯器進(jìn)行編譯。由于編譯器的差別，有少量語(yǔ)言規(guī)則不一致，所以需要作一些修改。如在gcc中，可以定義一個(gè)含有變量的可變長(zhǎng)度數(shù)組，但在MicrosoftVisualC++中則只能定義固定長(zhǎng)度的數(shù)組。當(dāng)然這僅是極小的差別，稍作修改以后程序是可移植的。本文使用渦格法源程序如下#include<stdio.h>//inputandoutputlibrary.#include<math.h>//compilewithgcc,youmustusetheargument"-lm"tolink<math.h>#definepi3.1415926#defineM300#defineN300intgetzp(doublep[],doublex[],doublez[],intnp,doublev0[],doubles,doublec);intinducedu(doubleppa[3],doubleppb[3],doubleppc[3],doubleiu[3]);intinduceduw(doubleppa[3],doublewvvector[3],doubleppc[3],doubleiuw[3]);intsolve(doublea[],doublesolu[],doubleb[],intn);intmain(){inti,j,k,l,I,J,m=15,n=20;double s=2.0,c=1.0;//s:spanlength;c:chordlength.doublev0[3],absv0=1.0;double vv0[3];//vv0:vectoroftheinletflowvelocity.// double a=0.0;// double a=2.5;// double a=5.0;double a=10.0;doublematrix[M*N],rhs[M],sol[M];doubleus[3],uc[3],uw[3],vinduced[3];//horseshoevortexinducedvelocity.double normal[3],pa[3],pb[3],pc[3];//normalvector,pointa,pointb,pointc.doublewvector[3];//wakevortexvector.//writetheresulttothefile!FILE*fp;/*//NACAa=0.8int ip,np=27;//ip:indexofpoint;np:numberofpoint.double xc[27]={0,0.005,0.0075,0.0125,0.025,0.05,0.075,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,0.975,1.0};double zf[27]={0,0.0423,0.0595,0.0907,0.1586,0.2712,0.3657,0.4482,0.5869,0.6993,0.7905,0.8635,0.9202,0.9615,0.9881,1.0,0.9971,0.9786,0.9434,0.8892,0.8121,0.7087,0.5425,0.3586,0.1713,0.0823,0};doublex[27],z[27],fmax;//x:x-coordinate,z:z-coordinate,fmax:maximumofthef.doubledx,dz,df;//fmax=0.02*c;fmax=0.04*c;for(ip=0;ip<np;ip++){x[ip]=c*xc[ip];z[ip]=fmax*zf[ip];}*//*//arcpropeller!int ip,np=51;//ip:indexofpoint; np:numberofpoint.doublex[51],z[51],fmax;//x:x-coordinate,z:z-coordinate,fmax:maximumofthef.double dx,dz,df;double x0,z0,r;//fmax=0.02*c;fmax=0.04*c;r=(c*c/4.0+fmax*fmax)/(2.0*fmax);x0=c/2.0;z0=fmax-r;for(ip=0;ip<np;ip++){x[ip]=c*ip/(np-1.0);z[ip]=z0+sqrt(r*r-(x[ip]-x0)*(x[ip]-x0));}*///parabolapropellerint ip,np=51;//ip:indexofpoint; np:numberofpoint.doublex[51],z[51],fmax;//x:x-coordinate,z:z-coordinate,fmax:maximumofthef.double dx,dz,df;double a0,a1,a2;//fmax=0.02*c;fmax=0.04*c;a0=0;a1=4.0*fmax/c;a2=-4.0*fmax/(c*c);for(ip=0;ip<np;ip++){x[ip]=c*ip/(np-1.0);z[ip]=a0+a1*x[ip]+a2*x[ip]*x[ip];}//Setvaluestotheinletflowvelocityvectorandthevectoroftheinflowvelocity.v0[0]=absv0*cos(a*pi/180);v0[1]=0.0;v0[2]=absv0*sin(a*pi/180);vv0[0]=cos(a*pi/180);vv0[1]=0.0;vv0[2]=sin(a*pi/180);//Setvaluestothematrix[][],rhs[]andsol[]beforeyouusethem,oryouwillgetweirdvalues.for(I=0;I<M;I++){for(J=0;J<N;J++){matrix[I*M+J]=0.0;}rhs[I]=0.0;sol[I]=0.0;}//Computetheelementvaluesofthematrixandtheright-hand-sidevector.//I:indicatetheinducedvelocitypoint;J:indicatethevortexsegment.for(I=0;I<M;I++){i=I/n;j=I-i*n;//i,j:indicatethemeshlocationoftheinducedvelocitypoint//computethecoordinateofthepointc.pc[0]=(j+0.75)*(c/n);pc[1]=(i+0.5)*(s/m);getzp(pc,x,z,np,vv0,s,c);//getthez-coordinateofthepoint.//Computetheright-hand-sidevector.for(ip=0;ip<np;ip++){if(pc[0]>=x[ip]&&pc[0]<x[ip+1]){dx=x[ip+1]-x[ip];dz=z[ip+1]-z[ip];df=sqrt(dx*dx+dz*dz);normal[0]=-dz/df;//note:itis"-dz/df",not"dz/df".normal[1]=0;normal[2]=dx/df;}}for(k=0;k<3;k++){rhs[I]=rhs[I]-v0[k]*normal[k];}for(J=0;J<N;J++){vinduced[0]=vinduced[1]=vinduced[2]=0.0;i=J/n;j=J-i*n;//i,j:indicatethemeshlocationofthevortexsegment.//computethevortexinducedvelocityiny-coordinate.//computethecoordinateofthepointa.pa[0]=(j+0.25)*(c/n);pa[1]=i*(s/m);getzp(pa,x,z,np,vv0,s,c);//getthez-coordinateofthepoint.//computethecoordinateofthepointb.pb[0]=(j+0.25)*(c/n);pb[1]=(i+1)*(s/m);getzp(pb,x,z,np,vv0,s,c);//getthez-coordinateofthepoint.inducedu(pa,pb,pc,us);//computeus.for(k=0;k<3;k++){vinduced[k]=vinduced[k]+us[k];}//computeuc(i+1).for(l=j;l<2*n;l++){//papa[0]=(l+0.25)*(c/n);pa[1]=(i+1)*(s/m);getzp(pa,x,z,np,vv0,s,c);//pbpb[0]=(l+1+0.25)*(c/n);pb[1]=(i+1)*(s/m);getzp(pb,x,z,np,vv0,s,c);inducedu(pa,pb,pc,uc);//computeuc(i+1).for(k=0;k<3;k++){vinduced[k]=vinduced[k]+uc[k];}//computeuc(i).for(l=j;l<2*n;l++){//papa[0]=(l+0.25)*(c/n);pa[1]=i*(s/m);getzp(pa,x,z,np,vv0,s,c);//pbpb[0]=(l+1+0.25)*(c/n);pb[1]=i*(s/m);getzp(pb,x,z,np,vv0,s,c);inducedu(pa,pb,pc,uc);//computeuc(i).for(k=0;k<3;k++){vinduced[k]=vinduced[k]-uc[k];}}//computeuw(i+1).//papa[0]=(2*n+0.25)*(c/n);pa[1]=(i+1)*(s/m);getzp(pa,x,z,np,vv0,s,c);//setthewakevortexvector.wvector[0]=vv0[0];wvector[1]=vv0[1];wvector[2]=vv0[2];induceduw(pa,wvector,pc,uw);//computeuw(i+1).for(k=0;k<3;k++){vinduced[k]=vinduced[k]+uw[k];}//computeuw(i).//papa[0]=(2*n+0.25)*(c/n);pa[1]=i*(s/m);getzp(pa,x,z,np,vv0,s,c);//wvectorwvector[0]=vv0[0];wvector[1]=vv0[1];wvector[2]=vv0[2];induceduw(pa,wvector,pc,uw);//computeuw(i).for(k=0;k<3;k++){vinduced[k]=vinduced[k]-uw[k];}//Setthematrixvalues.for(k=0;k<3;k++){matrix[I*M+J]=matrix[I*M+J]+vinduced[k]*normal[k];}}//matchesfor(J=0;J<N;J++){}//matchesfor(I=0;I<N;I++){printf("\ntheright-hand-sidevectoris:\n");for(I=0;I<M;I++){printf("(%d)%g",I,rhs[I]);}//SolvethelinearsystemwithSORiterationmethod.solve(matrix,sol,rhs,M);fp=fopen("4p100.txt","w");//themeaningofthename:'2'meansfmax=0.02*c;'n'meansNACAa=0.8propeller;'00'meanstheangleoftheinletvelocitya=0.0//inthesameway,'4'meansfmax=0.04*c;'a'meansarcpropeller;'p'meansparabolapropeller;//'25'meansa=2.5;'50'meansa=5.0;'100'meansa=10.0;//writethehorseshoevortexintensitydistributionofthemiddlepropelleriny-coordinate.for(I=((m-1)/2)*n;I<((m-1)/2+1)*n;I++){ //note:itis((m-1)/2)*n,not((m+1)/2)*n.fprintf(fp,"%f%g\n",(I%n+0.75)/n,sol[I]);}fclose(fp);return0;}//Ggetzp(doublep[],doublex[],doublez[],intnp,doublevv0[],doubles,doublec){if(p[0]>=0.0&&p[0]<=c){int ip;//ip:indexofpoint.double dx,dz;for(ip=0;ip<np;ip++){if(p[0]>=x[ip]&&p[0]<x[ip+1]){dx=x[ip+1]-x[ip];dz=z[ip+1]-z[ip];p[2]=z[ip]+dz*(p[0]-x[ip])/dx;//interpolatingthez-coordinateofthepoint.}}}else{inti;doublealpha;//avariablewhosevalueisbetween0and1.doublevt[3],vw[3];doublea[3],b[3];doubledx,dz,df;alpha=2*(p[1]/s-0.5)*(p[1]/s-0.5);//inordertodiscriptionthewakevortexdeformation.dx=x[np-1]-x[np-2];dz=z[np-1]-z[np-2];df=sqrt(dx*dx+dz*dz);vt[0]=dx/df;vt[1]=0.0;vt[2]=dz/df;for(i=0;i<3;i++){vw[i]=alpha*vv0[i]+(1-alpha)*vt[i];}b[0]=vw[2]/vw[0];//computetheslopeofthevector.b[1]=vv0[2]/vv0[0];b[2]=0.0;//solvethelinearsystembyanaliticalmethod.a[0]=(b[1]-b[0])/(2*c);a[1]=b[1]-4*a[0]*c;a[2]=b[2]-a[0]*c*c-a[1]*c;p[2]=a[0]*p[0]*p[0]+a[l]*p[0]+a[2];//z=al*xA2+a2*x+a3}return0;}//computethehorseshoevortexinducedvelocityonthesufaceofthepropeller,inducedu(doubleppa[3],doubleppb[3],doubleppc[3],doubleiu[3])//note:abisvortex,pointcisinducedvelocitypoint.{doublea,b,c,d,e,v;doubleavector[3],cvector[3],vvector[3];//vvector=(avectorXcvector)/(a*d).a=sqrt(pow(ppb[0]-ppa[0],2)+pow(ppb[l]-ppa[l],2)+pow(ppb[2]-ppa[2],2));//a=||pb-pa||b=sqrt(pow(ppb[0]-ppc[0],2)+pow(ppb[l]-ppc[l],2)+pow(ppb[2]-ppc[2],2));//b=||pb-pc||c=sqrt(pow(ppa[0]-ppc[0],2)+pow(ppa[l]-ppc[l],2)+pow(ppa[2]-ppc[2],2));//c=||pa-pc||e=(a*a+c*c-b*b)/(2*a);d=sqrt(c*c-e*e);//thedistancefromthepointctothevortexsegment.//Computetheabsolutevalueoftheinducedvelocity.//note:d=sqrt(c*c-e*e),maybed>=0.0l*cisareasonablejudgementcriteria!if(d>=0.0l*c){v=((e/c)+(a-e)/b)/(4*pi*d);}elseif(e<0){v=(l.0/(e*e)-l.0/pow(a-e,2))*d/(8*pi);//note:itis"l.0/(e*e)",not"l.0/e*e"!}else{v=(l.0/pow(a-e,2)-l.0/(e*e))*d/(8*pi);//note:itis"l.0/(e*e)",not"l.0/e*e"!}//Computethevectoroftheinducedvalue.//note:"avectorXcvector"isequalto"avectorXdvector".avector[0]=ppb[0]-ppa[0];avector[1]=ppb[1]-ppa[1];avector[2]=ppb[2]-ppa[2];//ppa-->ppbcvector[0]=ppc[0]-ppa[0];cvector[1]=ppc[1]-ppa[1];cvector[2]=ppc[2]-ppa[2];//ppa-->ppc//vvector=(avectorXcvector)/(a*d).note:"X"(multiplicationcross,crossproduct)vvector[0]=(avector[1]*cvector[2]-avector[2]*cvector[1])/(a*d);vvector[1]=-(avector[0]*cvector[2]-avector[2]*cvector[0])/(a*d);vvector[2]=(avector[0]*cvector[1]-avector[1]*cvector[0])/(a*d);iu[0]=v*vvector[0];iu[1]=v*vvector[1];iu[2]=v*vvector[2];return0;}//Subroutineinduceduw():induceduw(doubleppa[3],doublewvvector[3],doubleppc[3],doubleiuw[3]){intk;doublec,d,e,v;doublecvector[3],vvector[3];c=sqrt(pow(ppa[0]-ppc[0],2)+pow(ppa[1]-ppc[1],2)+pow(ppa[2]-ppc[2],2));//c=||pa-pc||cvector[0]=ppc[0]-ppa[0];cvector[1]=ppc[1]-ppa[1];cvector[2]=ppc[2]-ppa[2];//ppa-->ppce=0.0;for(k=0;k<3;k++){e+=cvector[k]*wvvector[k];//e=cvector[]multiplicationdot(innerproduct)wvvector[].}d=sqrt(c*c-e*e);//Computetheabsolutevalueoftheinducedvelocity.if(d>=0.01*c){v=(e/c+1)/(4*pi*d);}elseif(e<0){v=d/(8*pi*e*e);}else{v=(2-d*d/(2*e*e))/(4*pi*d);}//Computethevectoroftheinducedvalue.//vvector=(wvvectorXcvector)/d.note:"X"(multiplicationcross,crossproduct)vvector[0]=(wvvector[1]*cvector[2]-wvvector[2]*cvector[1])/d;vvector[1]=-(wvvector[0]*cvector[2]-wvvector[2]*cvector[0])/d;vvector[2]=(wvvector[0]*cvector[1]-wvvector[1]*cvector[0])/d;//Computethewakevortexinducedvelocity.iuw[0]=v*vvector[0];iuw[1]=v*vvector[1];iuw[2]=v*vvector[2];return0;}//SolvethelinearsystemwithSORsolve(doublea[],doublesolu[],doubleb[],intn){inti,j,k,NN=1000;//now"N"isaconstant.doublee=1e-5,w=1.0;//e:thecriteriaoftheconvergence;w:relaxationvariable.doublet[M];//cannotallocateanarrayofconstantsize0 >t[n]doublesum=0.0;//m=0;for(k=0;k<NN&&m!=n;k++){m=0