理想流體力學大作業(yè)學生姓名：學號：2023年10月Hess—Smith方法計算物體附加質量摘要：本文運用Hess-Smith方法計算了圓球、橢球和圓柱的附加質量系數(shù)以及橢球并行的干擾效應。同時，文章分析了網(wǎng)格變化對計算值的影響趨勢。本文使用matlab語言對圓球、橢球與圓柱的模型進行了網(wǎng)格有限元的劃分，得到各個單元的節(jié)點坐標，然后利用Hess-Smith方法對圓球、橢球及并行橢球的附加質量系數(shù)進行計算及分析。關鍵字：邊界元；Hess-Smith;附加質量系數(shù)一、物理背景Hess-Smith方法是一種計算任意三維物體勢流的方法，該方法由美國的Hess和Smith兩人于20世紀60年代提出。Hess-Smith方法又稱為分布奇點法，作為一種邊界元方法，它用許多平面四邊形或三角形外表單元來表示物體外表，并在每個單元上布置強度未知的源，然后在物體外表的某些考察點上滿足法向速度為零的物面邊界條件，得到求單元源密度的線性代數(shù)方程組。求解方程組得到源密度分布，進而可求流場內任意點的速度、壓力等物理量。二、理論依據(jù)2.1分布源模型的建立為無界流中的物體外表，來流為均勻流，在無窮遠處流體的速度為：()為定常速度勢，并在物體外部空間域中滿足拉普拉斯方程，在物面上適合不可穿透條件，在無窮遠處，應該與均勻來流的速度勢相同。即〔物體外〕()〔物面上〕()其中，單位法線向量指向物體內部。在速度勢中分出的均勻來流項，記()這里的是擾動速度勢，應適合以下定解條件：()用Rpq表示點p和點q之間的距離，根據(jù)格林第三公式，當p點位于物面s外部和遠方控制面c的內部之空間域時，有如下公式：()由遠方邊界條件可知，遠方封閉控制面上的積分趨于零，從而上式化為：()又由式〔〕可得：()得到混合分布模型，為了得到單一分布模型表示的擾動勢，在物體內部域中構造一個適宜的內部解。于上述物體外部的點P，函數(shù)在物體內部域中沒有奇點，在物體內部域中對函數(shù)和用格林第三公式，得：()將〔〕與〔〕相減，得：()取下式定解條件中的：(在內部)()(在)()那么式()成為：()其中，()2.2分布源密度的求解式()中右端分布源的法向導數(shù)極限由兩局部組成，一局部是p點附近小曲面ε的奉獻，另一局部是物面其余局部的奉獻。法向指向取向物體內部，小曲面ε的奉獻為2πσ〔p〕,那么有如下關系式：()再結合物面條件()，得到()這就是分布源密度所適合的線性積分方程。把積分方程()轉換成線性代數(shù)方程組，即用離散量代替連續(xù)變量。把物面分成小塊，記()用平面四邊形或三角形來近似代替小曲面。具體做法如下，取第小塊的四個頂點坐標之算術平均值，得到中心點的坐標。計算對角線向量的向量積〔指向與曲面法線指向相符合〕，用表示該方向上的單位向量，形成以為法線且通過中心點的平面，再把四個頂點向該平面作投影，以四個投影點為頂點組成平面四邊形，用代替原來的小曲面，稱為單元。通常把小范圍內的分布源密度作為常數(shù)，因此只要分割不太粗，可以認為在單元上為常數(shù)，記作，從而()因此物面上的積分可以用個平面四邊形〔三角形〕上積分之和來近似，即()上式左端的未知量是連續(xù)型變量，而上式右端的未知量是個離散量。為了求解這個未知數(shù)，須要個方程。取積分方程()中的動點為個單元的中心點，稱之為控制點，即控制物面條件使之成立的點。用近似式(2.2.5)代替積分方程(2.2.2)的左端，便可以寫出的階線性代數(shù)方程組：()當計算出影響系數(shù)后，即可解線性方程組得到分布源密度。2.3速度勢與附加質量的求解根據(jù)速度勢在控制點pi處的值，由公式：()()根據(jù)2.2得到的分布源密度，求解線性方程組〔〕可得速度勢的值。物體的附加質量，表示物體沿方向運動引起的方向的附加質量，公式如下:()()根據(jù)所求得的速度勢的值可計算處附加質量的值。三、數(shù)值模型及參數(shù)計算3.1數(shù)值模型要求解流場中物體外表的速度勢分布，需要先將物體的外外表進行網(wǎng)格劃分。經(jīng)過網(wǎng)格劃分以后，原來的物體連續(xù)外外表被離散為N×M個相對獨立的小平面，這些小平面構成了求解該問題的數(shù)值模型。Hess-Smith的根本思想是將連續(xù)曲面的積別離散為小單元來簡化計算，其計算思路核心在于解該方程組：，通過求解線性方程得到。對于不同的計算目的，只需要改變控制面條件，即改變來實現(xiàn)。得到后，進而由求及附加質量，其中：為求，令，那么求得。故而，同理可得。3.2參數(shù)計算3.2.1計算對于球面：對于橢球面：設，，那么易知，，，故?！瞐=1即為球〕對于圓柱面：在圓柱頂面上，同理在其底面上有，側面上那么有3.2.2計算對于上述幾何模型中的四邊形上均勻分布奇點的誘導速度公式計算，首先將四邊形上的分布源密度取為1，四邊形四個頂點逆時針方向排列，頂點pi的坐標為〔ξi，ηi，0〕，其中坐標系oξηζ的原點為四邊形形心，平面oξη即為四邊形平面，那么需計算如下三個積分式：各積分式的計算公式為：在本幾何模型中，需要將世界坐標系轉換到局部坐標系中，查閱相關文獻可知，當給定物體上的3個不共線的點，即可建立局部坐標系。其中局部坐標系的原點為點，軸的正向為的方向，軸的正向為的方向，軸的正向為的方向。設軸的單位矢量分別為它們在世界坐標系中的方向余弦分別用表示。對軸，可按下式給出：，，，為世界坐標系3根軸的單位矢量。矩陣就是由世界坐標系到物坐標系的坐標變換矩陣。經(jīng)過坐標轉換可得每個網(wǎng)格單元的Sx、Sy、Sz，然后得到世界坐標系下的S，即由于編程實現(xiàn)上述過程非常繁復，現(xiàn)采用一種更為簡便方法的求解當時，當時，，從而求得的系數(shù)矩陣。3.2.3計算由可得當時，當時，如圖，設位于上，取一小圓ε，轉換為極坐標積分：〔廣義積分收斂〕故可近似，得到系數(shù)矩陣。3.2.4計算m11、m33令，那么同理求得。對于雙橢球體，如下列圖所示，只需將第二個橢球的單元標記為〔，…〕即可。方程組的形式不變，在計算m11或m33時分別對兩個橢球面的單元進行積分即可分別得到m11a、m11b、m33a、m33b。四、幾何模型對于連續(xù)的積分方程，通常的數(shù)值處理是轉換成線性代數(shù)方程組，即用離散量代替連續(xù)變量。在幾何建模上，就是將物面s分成N小塊，記為采用平面四邊形或三角形近似代替小曲面。4.1球面網(wǎng)格劃分4.1.1網(wǎng)格劃分使用球坐標對球體進行網(wǎng)格劃分，取球體半徑R=1，以球心O為原點建立坐標系如下列圖所示取xoz平面上的圓，圓上任意一點與x軸的夾角范圍為0-π，劃分N份，每相鄰兩份夾角為δφ=π/N，在球體上以球心為圓心，劃分N條緯線。緯線為在球體上的同心圓，取任意同心圓，圓半徑為r，將圓劃分M份，每相鄰兩份之間夾角為δθ=2π/M,N個同心圓都劃分M份，形成M條經(jīng)線。由縱橫交錯的線將球體劃分形成網(wǎng)格。取球體任意網(wǎng)格交點A,設坐標為A(x，y，z),由球坐標可將A坐標表示為：在matlab中根據(jù)以上網(wǎng)格劃分原理劃分網(wǎng)格，得到的球體幾何模型為：4.1.2節(jié)點坐標取球體面上任意一個小網(wǎng)格，在小網(wǎng)格ABCD上，A標記為A(i，j)，即縱向第i個點，橫向第j個點，那么i=1,2,…..,N+1j=1,2,……,M+1B,C,D可由A(i，j)表示，如上圖所示。即可得到個小網(wǎng)格單元，故有個控制點P。對于頂點處，單元為三角形：或是4.1.3球面網(wǎng)格小單元面積對于四邊形單元，AB//CD,由對稱性可知，ABCD為等腰梯形。故ABCD為平面。ABCD的面積為δS。圖解如圖：梯形面積公式為：其中梯形高近似為腰長計算，得如下公式：EF的長度計算公式如下：梯形面積可表示為：4.2橢球面網(wǎng)格劃分4.2.1網(wǎng)格劃分同球形網(wǎng)面劃分類似，得到N×M個小網(wǎng)格單元。與球面不同的是，沿x軸方向延長。由matlab劃分網(wǎng)格得到的橢球幾何圖形如下列圖所示。4.2.2節(jié)點坐標長短軸比為d。橢球面上任意網(wǎng)格交點坐標為:那么橢球體外表每個單元的控制點坐標為：4.2.3網(wǎng)格單元面積在O-XYZ坐標系中，橢球體的剖面圖如下列圖所示，易知：橢球面可由球面沿x軸延伸得到，即橢球面上某一小網(wǎng)格面ABCD的AC邊延長為l2，由球面小網(wǎng)格面可知，面積為原來的l2/l1。即：即得小網(wǎng)格面積。4.2.4面積公式驗證由3.2.3節(jié)得到了每個網(wǎng)格單元面積，由此可以得到橢球體外表積的近似值，為了減少計算過程中的誤差，查閱相關文獻后得到橢球體外表積計算公式，然后驗證兩者誤差。設三橢球面方程的方程為：。其面積為：這里的，，又，分別為第一種和第二種橢球積分。對于，，k=0，那么F=u，E=u。當d=5，N=40，M=32時，計算出的單元面積和，理論值。相對誤差為，誤差很小，說明公式的正確性與精確性。4.3圓柱面網(wǎng)格劃分4.3.1網(wǎng)格劃分使用柱坐標對圓柱體進行網(wǎng)格劃分，取圓柱底面半徑R=1，以底面圓心O為原點建立坐標系。取yoz平面上的圓，圓上任意一點與y軸的所稱角度范圍為0-2π，劃分K份，每相鄰兩份夾角為δθ=2π/N，在圓柱上按等距將圓柱沿高度方向分為M段，每份高δL那么由縱橫交錯的線可將圓柱側外表劃分形成網(wǎng)格。從而圓柱側面可分為M*K個面元。按圓柱上下底面的極徑均分為N份，每份長δr那么同理可以將兩個底面分為2N*K個面元，總計將圓柱分為(2N+M)*K個面元。在matlab中根據(jù)以上網(wǎng)格劃分原理劃分網(wǎng)格，得到的柱體幾何模型為：4.3.2節(jié)點坐標取圓柱任意網(wǎng)格交點A,設坐標為A(x，y，z),由柱坐標可將A坐標表示為：那么橢球體外表每個單元的控制點坐標為：4.3.3圓柱面網(wǎng)格小單元面積對于圓柱側外表，AB//CD,AC//BD,易知ABCD為矩形。其中矩形沿圓柱高度方向邊長為a=δL，另一邊邊長那么為對于圓柱頂面和底面可以按照兩三角形面積之差計算五、計算結果顯示5.1圓球附加質量計算結果經(jīng)過不同組〔改變網(wǎng)格數(shù)量〕數(shù)值模型的求解，得到圓球附加質量系數(shù)M11、M33的計算結果，如下表4.1所示。表4.1單球體附加質量系數(shù)隨球面網(wǎng)格數(shù)量的變化情況網(wǎng)格數(shù)50200450648968M110.60610.54550.52670.52110.5163M330.56290.52540.51620.51330.5108為了更加直觀地考察附加質量隨球面網(wǎng)格數(shù)量的變化情況，繪制出與上表4.1相應的折線圖，這里將M11和M33的變化曲線畫在同一張圖里，便于觀察比照，如下列圖4.1所示。圖4.1單個球體附加質量系數(shù)隨網(wǎng)格數(shù)變化情況由上圖4.1顯示的曲線變化趨勢可以看出，隨著網(wǎng)格數(shù)的增加，數(shù)值計算的結果M11和M33都收斂于理論的解析解0.5。另外，所有的數(shù)值計算結果都大于理論值0.5。本項的計算結果，驗證了球體外表在流場中的各向同性的性質，即對于球體，有M11=M33，且在一定程度上驗證了解析解為0.5的正確性?；蛘叻催^來講，數(shù)值計算的結果收斂于解析解，反映出HESS-SMITH方法的可行性和正確性。下面，從另一個角度來考察附加質量的性質。即改變物體模型，考察附加質量在各個方向上表現(xiàn)出來的性質如何。5.2橢球附加質量計算結果經(jīng)過不同組〔改變網(wǎng)格數(shù)量〕數(shù)值模型的求解，得到橢球附加質量系數(shù)M11、M33的計算結果，如下表4.2所示。表4.2橢球附加質量系數(shù)隨網(wǎng)格數(shù)量的變化情況網(wǎng)格數(shù)320490640980M110.06310.06410.06310.0608M330.92990.91540.93010.8942繪制出與上表4.2相應的折線圖。這里因為M11和M33在量級上相差較大，故將M11和M33的變化曲線畫在兩張不同的圖里，以便于觀察各自的走向和變化。M11和M33的變化曲線分別如下列圖4.2、圖4.3所示。圖4.2橢球附加質量系數(shù)M11隨網(wǎng)格數(shù)變化情況圖4.3橢球附加質量系數(shù)M33隨網(wǎng)格數(shù)變化情況5.30雙橢球附加質量計算結果對于計算兩個橢球的干擾效應，計算時每個橢球劃分為980個單元，軸間距按照短軸的3、5、7倍分別計算,得到兩個橢球的附加質量系數(shù)。計算結果如表4.3所示：表4.3雙橢球橢球附加質量系數(shù)隨網(wǎng)格數(shù)量的變化情況間距比357M11a0.07610.06770.0644M11b0.07610.06680.0629M33a0.94890.86710.8493M33b0.94590.86300.8470依據(jù)表4.3繪制出圖4.4及4.5圖4.4雙橢球M11隨軸間距變化關系圖4.5雙橢球M33隨軸間距變化關系5.4圓柱附加質量計算結果經(jīng)過不同組〔改變網(wǎng)格數(shù)量〕數(shù)值模型的求解，得到圓球附加質量系數(shù)M11、M33的計算結果，如下表4.4所示。表4.4圓柱附加質量系數(shù)隨球面網(wǎng)格數(shù)量的變化情況網(wǎng)格數(shù)375525750900M110.13420.13040.13460.1352M330.89820.90250.90500.9148繪制出與上表4.4相應的折線圖。M11和M33的變化曲線分別如下列圖4.6、圖4.7所示。圖4.6圓柱附加質量系數(shù)M11隨網(wǎng)格數(shù)變化情況圖4.7圓柱附加質量系數(shù)M33隨網(wǎng)格數(shù)變化情況六、結果討論誤差分析6.1結果討論對于單個球體，在無界流中其附加質量系數(shù)理論值是0.5，從計算結果及圖表能看出，當球面單元劃分數(shù)增加時，M11、M33均收斂于0.5。但M11、M33均大于0.5，這說明即便單元劃分數(shù)目較大，但其外表仍不是光滑曲面，相比于球面，由水動力的知識易于理解其附加質量系數(shù)有微幅增量。圖表知，當球面面元數(shù)為968時，M11、M33相對誤差分別約為3.26%、2.16%，誤差在允許范圍內，初步驗證了自編程序的正確性可可靠性。對于單個橢球，可從蘭伯附加質量系數(shù)圖表查得，當d=5時，其M11、M33理論值分別約為0.053、0.88,取面元數(shù)為980時的結果進行比擬，相對誤差分別為14.72%、1.61%。M33的誤差已經(jīng)很小，可以認為計算結果相當準確。對于M11相對誤差較大的原因稍后論述。對于雙橢球，從圖表能直觀的看出兩個橢球之間的擾動效應。當橢球附近流域內有其他物體時，其附加質量系數(shù)會增大，增量隨間距為遞減。理論上當兩個橢球間距無窮遠時，單個橢球的附加質量即為無界流中的附加質量，但從圖中可以看出當軸間距為7倍以上時，單個橢球的附加質量就非常接近無擾動的值，此時可初步認為7倍間距為臨界距離。另外理論分析知，兩個橢球的附加質量系數(shù)應相等，從圖4.4、圖4.5可以直觀了解到兩個橢球的附加質量系數(shù)誤差很小，與理論結果相符。對于圓柱，其M11的值未知，M33的理論值為1。從計算結果可以看出，隨著網(wǎng)格劃分數(shù)目的增加，M33的值逐漸趨近于1，當網(wǎng)格劃分至900個面元時M33的相對誤差為8.52%，與真實值比擬接近。推測其誤差較大是因為網(wǎng)格劃分過程中網(wǎng)格大小不均勻所致，位于圓柱上下底面圓心附近處的面元過于細長影響了運算結果的準確性。6.2誤差原因分析1.方法誤差〔截斷誤差〕在模型建立中，以小平面代曲面，本身就帶來誤差。求δS時，GH,EF的均舍去了高階項O(x2),帶來的極小誤差。求Cij時，對于i=j,Cij≈2〔πδS〕1/2，由于是近似解，對于φ〔P〕的精準性有一定影響。此外模型建立過程中，網(wǎng)格劃分如果長寬比過大，會使Cij的誤差更為嚴重，而對于橢球而言，由于其結構較為細長，因此會產生大量細長的三角形單元；對于圓柱而言，在其底面上也會因網(wǎng)格劃分過密導致這一問題，使得方法誤差的影響加大。2.舍入誤差在計算過程中由于計算機存儲位數(shù)固定，數(shù)值在計算中會不停地在規(guī)定位數(shù)上取舍，帶來舍入誤差。即便計算時數(shù)值取雙精度，但由于參與計算的參數(shù)數(shù)量大，計算過程長，造成誤差積累，而求解過程中由于涉及到線性方程組的求解，因而舍入誤差的量不容無視。由于M11量級相對較小，此時舍入誤差占總誤差的比值增大，因而此時會導致M11相對誤差到達12%。6.3減小誤差措施改變網(wǎng)格結構，使單元盡量接近正方形；適當增加網(wǎng)格數(shù)，增加方法近似性，減小截斷誤差。但數(shù)目不宜過大，會導致程序占用內存過大，計算速度明顯下降，舍入誤差累和增大等不利影響。參考文獻[1]戴遺山,段文洋.船舶在波浪中運動的勢流理論.國防工業(yè)出版社,2006.[2]范尚雍.船舶操縱性.國防工業(yè)出版社.1988[3]張亮，李云波.流體力學.哈爾濱工程大學出版社.2006.附錄：matlab自編程序[1]圓〔橢〕球附加質量計算[2]圓柱附加質量計算[3]雙橢球附加質量干擾計算[1]圓〔橢〕球附加質量計算N=input('輸入N：');M=input('輸入M：');R=1;%單元劃分數(shù)N*M,球半徑為1；d=input('輸入d：');%長短軸比%劃分結構化網(wǎng)格phi=linspace(0,pi,N+1);thet=linspace(0,2*pi,M+1);alpha=zeros(1,N);x=zeros(N+1,M+1);y=zeros(N+1,M+1);z=zeros(N+1,M+1);P=cell(N,M);A=cell(N+1,M+1);delta_S=zeros(N,M);delta_phi=pi/N;delta_thet=2*pi/M;fori=1:Nalpha(i)=acos(((sin(i*delta_phi)-sin((i-1)*delta_phi)))/delta_phi);endfori=1:N+1;forj=1:M+1;y(i,j)=R*sin(phi(i))*cos(thet(j));z(i,j)=R*sin(phi(i))*sin(thet(j));x(i,j)=d*R*cos(phi(i));A{i,j}=[x(i,j),y(i,j),z(i,j)];endendforj=1:M;P{1,j}=1/3*(A{1,j}+A{2,j}+A{2,j+1});P{N,j}=1/3*(A{N,j}+A{N,j+1}+A{N+1,j});enddelta_S(1,:)=cos(alpha(1))/cos(atan(d*tan(alpha(1))))*R^2*delta_phi*delta_thet*sin(delta_phi*(0.5));delta_S(N,:)=delta_S(1,1);fori=2:N-1;delta_S(i,:)=abs(cos(alpha(i))/cos(atan(d*tan(alpha(i)))))*R^2*delta_phi*delta_thet*sin(delta_phi*(i-0.5));endfori=2:N-1;forj=1:M;P{i,j}=0.25*(A{i,j}+A{i,j+1}+A{i+1,j}+A{i+1,j+1});endendsurf(x,y,z,x)xlabel('x')ylabel('y')zlabel('z')axisimage%開始進行單元計算，解方程組%設遠方速度為1，x正方向a=zeros(N*M,N*M);b=zeros(N*M,1);c=zeros(N*M,N*M);sigma=zeros(N*M,1);PHI=zeros(N*M,1);r=zeros(N*M,N*M);%求橢球m_11n_p=cell(N*M,1);PHI_1=zeros(N*M,1);fork=1:N*Mn_p{k}(1)=-1/d^2*P{k}(1)/sqrt(P{k}(1)^2/d^4+P{k}(2)^2+P{k}(3)^2);b(k)=-1*n_p{k}(1);forl=1:N*Mr(k,l)=sqrt((P{k}(1)-P{l}(1))^2+(P{k}(2)-P{l}(2))^2+(P{k}(3)-P{l}(3))^2);a(k,l)=delta_S(l)*((P{k}(1)-P{l}(1))*P{k}(1)/d^2+(P{k}(2)-P{l}(2))*P{k}(2)+(P{k}(3)-P{l}(3))*P{k}(3))/r(k,l)^3/sqrt(P{k}(1)^2/d^4+P{k}(2)^2+P{k}(3)^2);c(k,l)=delta_S(l)/r(k,l);enda(k,k)=2*pi;c(k,k)=2*sqrt(pi*delta_S(k));endsigma=linsolve(a,b);m_11=0;fork=1:N*Mforl=1:N*MPHI_1(k)=PHI_1(k)+c(k,l)*sigma(l);endendfork=1:N*Mm_11=m_11-n_p{k}(1)*(PHI_1(k))*delta_S(k);end%求橢球m_33PHI_3=zeros(N*M,1);fork=1:N*Mn_p{k}(3)=-P{k}(3)/sqrt(P{k}(1)^2/d^4+P{k}(2)^2+P{k}(3)^2);b(k)=-1*n_p{k}(3);endsigma=linsolve(a,b);m_33=0;fork=1:N*Mforl=1:N*MPHI_3(k)=PHI_3(k)+c(k,l)*sigma(l);endendfork=1:N*Mm_33=m_33-n_p{k}(3)*(PHI_3(k))*delta_S(k);endM11=m_11/(4/3*pi*d)M33=m_33/(4/3*pi*d)[2]圓柱附加質量計算N=input('輸入N：');M=input('輸入M：');K=input('輸入K：');R=1;%單元劃分數(shù)(2N+M)*K,圓柱底面半徑為1；%劃分結構化網(wǎng)格L=10*R;rr=linspace(0,R,N+1);ll=linspace(0,L,M+1);thet=linspace(0,2*pi,K+1);x=zeros(2*N+M+1,K+1);y=zeros(2*N+M+1,K+1);z=zeros(2*N+M+1,K+1);P=cell(2*N+M,K);A=cell(2*N+M+1,K+1);delta_S=zeros(2*N+M,K);delta_r=R/N;delta_l=L/M;delta_thet=2*pi/K;fori=1:N;x(i,:)=0;forj=1:K+1;z(i,j)=rr(i)*sin(thet(j));y(i,j)=rr(i)*cos(thet(j));A{i,j}=[x(i,j),y(i,j),z(i,j)];endendfori=N+M+1:2*N+M+1;x(i,:)=L;forj=1:K+1;z(i,j)=rr(2*N+M+2-i)*sin(thet(j));y(i,j)=rr(2*N+M+2-i)*cos(thet(j));A{i,j}=[x(i,j),y(i,j),z(i,j)];endendfori=N+1:N+M;x(i,:)=ll(i-N);forj=1:K+1;z(i,j)=R*sin(thet(j));y(i,j)=R*cos(thet(j));A{i,j}=[x(i,j),y(i,j),z(i,j)];endendforj=1:K;P{1,j}=1/3*(A{1,j}+A{2,j}+A{2,j+1});P{2*N+M,j}=1/3*(A{2*N+M,j}+A{2*N+M,j+1}+A{2*N+M+1,j});endfori=1:N;delta_S(i,:)=sin(delta_thet)*(rr(i+1)^2-rr(i)^2)/2;delta_S(2*N+M-i+1,:)=delta_S(i,:);endsw=sin(delta_thet/2)*R*2*delta_l;fori=N+1:N+M;delta_S(i,:)=sw;endfori=2:2*N+M-1;forj=1:K;P{i,j}=0.25*(A{i,j}+A{i,j+1}+A{i+1,j}+A{i+1,j+1});endendsurf(x,y,z,x)xlabel('x')ylabel('y')zlabel('z')axisimagePP=2*N+M;%開始進行單元計算，解方程組%設遠方速度為1，x正方向a=zeros(PP*K,PP*K);b=zeros(PP*K,1);c=zeros(PP*K,PP*K);sigma=zeros(PP*K,1);PHI=zeros(PP*K,1);r=zeros(PP*K,PP*K);%求圓柱m_11n_p=cell(PP*K,1);PHI_1=zeros(PP*K,1);fork=1:K;forj=1:N;t=(k-1)*PP+j;n_p{t}(1)=1;n_p{t}(2)=0;n_p{t}(3)=0;endny=-cos(thet(k)+pi/K);nz=-sin(thet(k)+pi/K);forj=N+1:N+M;t=(k-1)*PP+j;n_p{t}(1)=0;n_p{t}(2)=ny;n_p{t}(3)=nz;endforj=N+M+1:N*2+M;t=(k-1)*PP+j;n_p{t}(1)=-1;n_p{t}(2)=0;n_p{t}(3)=0;endendfork=1:PP*Kb(k)=-1*n_p{k}(1);forl=1:PP*Kr(k,l)=sqrt((P{k}(1)-P{l}(1))^2+(P{k}(2)-P{l}(2))^2+(P{k}(3)-P{l}(3))^2);a(k,l)=-delta_S(l)*((P{k}(1)-P{l}(1))*n_p{k}(1)+(P{k}(2)-P{l}(2))*n_p{k}(2)+(P{k}(3)-P{l}(3))*n_p{k}(3))/r(k,l)^3;c(k,l)=delta_S(l)/r(k,l);enda(k,k)=2*pi;c(k,k)=2*sqrt(pi*delta_S(k));endsigma=linsolve(a,b);m_11=0;fork=1:PP*Kforl=1:PP*KPHI_1(k)=PHI_1(k)+c(k,l)*sigma(l);endendfork=1:PP*Km_11=m_11-n_p{k}(1)*(PHI_1(k))*delta_S(k);end%求圓柱m_33PHI_3=zeros(PP*K,1);fork=1:PP*Kb(k)=-1*n_p{k}(3);forl=1:PP*Kr(k,l)=sqrt((P{k}(1)-P{l}(1))^2+(P{k}(2)-P{l}(2))^2+(P{k}(3)-P{l}(3))^2);a(k,l)=-delta_S(l)*((P{k}(1)-P{l}(1))*n_p{k}(1)+(P{k}(2)-P{l}(2))*n_p{k}(2)+(P{k}(3)-P{l}(3))*n_p{k}(3))/r(k,l)^3;c(k,l)=delta_S(l)/r(k,l);enda(k,k)=2*pi;c(k,k)=2*sqrt(pi*delta_S(k));endsigma=linsolve(a,b);m_33=0;fork=1:PP*Kforl=1:PP*KPHI_3(k)=PHI_3(k)+c(k,l)*sigma(l);endendfork=1:PP*Km_33=m_33-n_p{k}(3)*(PHI_3(k))*delta_S(k);endM11=m_11/pi/L/R^2M33=m_33/pi/L/R^2[3]雙橢球附加質量干擾計算N=input('輸入N：');M=input('輸入M：');R=1;%單元劃分數(shù)N*M,球半徑為1；d=5;%長短軸比e=input('輸入e：');%軸間距與短軸之比%劃分結構化網(wǎng)格phi=linspace(0,pi,N+1);thet=linspace(0,2*pi,M+1);alpha=zeros(1,N);x=zeros(N+1,2*(M+1));y=zeros(N+1,2*(M+1));z=zeros(N+1,2*(M+1));P=cell(N,2*M);A=cell(N,2*M);delta_S=zeros(N,2*M);delta_phi=pi/N;delta_thet=2*pi/M;fori=1:Nalpha(i)=acos(((sin(i*delta_phi)-sin((i-1)*delta_phi)))/delta_phi);endfori=1:N+1;forj=1:M+1;y(i,j)=R*sin(phi(i))*cos(thet(j));z(i,j)=R*sin(phi(i))*sin(thet(j));x(i,j)=d*R*cos(phi(i));A{i,j}=[x(i,j),y(i,j),z(i,j)];endforj=M+2:2*(M+1);y(i,j)=e*R+R*sin(phi(i))*cos(thet(j-(M+1)));z(i,j)=R*sin(phi(i))*sin(thet(j-(M+1)));x(i,j)=d*R*cos(phi(i));A{i,j}=[x(i,j),y(i,j),z(i,j)];endendforj=1:2*M;P{1,j}=1/3*(A{1,j}+A{2,j}+A{2,j+1});P{N,j}=1/3*(A{N,j}+A{N,j+1}+A{N+1,j});enddelta_S(1,:)=cos(alpha(1))/cos(atan(d*tan(alpha(1))))*R^2*delta_phi*delta_thet*sin(delta_phi*(0.5));delta_S(N,:)=delta_S(1,1);fori=2:N-1;delta_S(i,:)=abs(cos(alpha(i))/cos(atan(d*tan(alpha(i)))))*R^2*delta_phi*delta_thet*sin(delta_phi*(i-0.5));endfori=2:N-1;forj=1:2*M;P{i,j}=0.25*(A{i,j}+A{i,j+1}+A{i+1,j}+A{i+1,j+1});endendholdonsurf(x(:,1:M+1),y(:,1:M+1),z(:,1:M+1),x(:,1:M+1));surf(x(:,M+2:2*(M+1)),y(:,M+2:2*(M+1)),z(:,M+2:2*(M+1)),x(:,1:M+1));xlabel('x')ylabel('y')zlabel('z')axisimage%開始進行單元計算，解方程組%設遠方速度為1，x正方向a=zeros(N*2*M,N*2*M);b=zeros(N*2*M,1);c=zeros(N*2*M,N*2*M);sigma=zeros(N*2*M,1);PHI=zeros(N*2*M,1);r=zeros(N*2*M,N*2*M);%求橢球m_11n_p=cell(N*2*M,1);PHI_1=zeros(N*2*M,1);fork=1:N*Mn_p{k}(1)=-1/d^2*P{k}(1)/sqrt(P{k}(1)^2/d^4+P{k}(2)^2+P{k}(3)^2);b(k)=-1*n_p{k}(1);forl=1:2*N*Mr(k,l)=sqrt((P{k}(1)-P{l}(1))^2+(P{k}(2)-P{l}(2))^2+(P{k}(3)-P{l}(3))^2);a(k,l)=delta_S(l)*((P{k}(1)-P{l}(1))*P{k}(1)/d^2+(P{k}(2)-P{l}(2))*P{k}(2)+(P{k}(3)-P{l}(3))*P{k}(3))/r(k,l)^3/sqrt(P{k}(1)^2/d^4+P{k}(2)^2+P{k}(3)^2);c(k,l)=delta_S(l)/r(k,l);enda(k,k)=2*pi;c(k,k)=2*sqrt(pi*delta_S(k));endfork=N*M+1:2*N*Mn_p{k}(1)=n_p{k-N*M}(1);b(k)=-1*n_p{k}(1);forl=1:2*N*Mr(k,l)=sqrt((P{k}(1)-P{l}(1))^2+(P{k}(2)-P{l}(2))^2+(P{k}(3)-P{l}(3))^2);a(k,l)=delta_S(l)*((P{k}(1)-P{l}(1))*P{k}(1)/d^2+(P{k}(2)-P{l}(2))*(P{k}(2)-e*R)+(P{k}(3)-P{l}(3))*P{k}(3))/r(k,l)^3/sqrt(P{k}(1)^2/d^4+(P{k}(2)-e*R)^2+