哈尔滨工程大学-理想流体力学-大作业

理想流体力学大作业学生姓名：学号：2023年10月Hess—Smith方法计算物体附加质量摘要：本文运用Hess-Smith方法计算了圆球、椭球和圆柱的附加质量系数以及椭球并行的干扰效应。同时，文章分析了网格变化对计算值的影响趋势。本文使用matlab语言对圆球、椭球与圆柱的模型进行了网格有限元的划分，得到各个单元的节点坐标，然后利用Hess-Smith方法对圆球、椭球及并行椭球的附加质量系数进行计算及分析。关键字：边界元；Hess-Smith;附加质量系数一、物理背景Hess-Smith方法是一种计算任意三维物体势流的方法，该方法由美国的Hess和Smith两人于20世纪60年代提出。Hess-Smith方法又称为分布奇点法，作为一种边界元方法，它用许多平面四边形或三角形外表单元来表示物体外表，并在每个单元上布置强度未知的源，然后在物体外表的某些考察点上满足法向速度为零的物面边界条件，得到求单元源密度的线性代数方程组。求解方程组得到源密度分布，进而可求流场内任意点的速度、压力等物理量。二、理论依据2.1分布源模型的建立为无界流中的物体外表，来流为均匀流，在无穷远处流体的速度为：()为定常速度势，并在物体外部空间域中满足拉普拉斯方程，在物面上适合不可穿透条件，在无穷远处，应该与均匀来流的速度势相同。即〔物体外〕()〔物面上〕()其中，单位法线向量指向物体内部。在速度势中分出的均匀来流项，记()这里的是扰动速度势，应适合以下定解条件：()用Rpq表示点p和点q之间的距离，根据格林第三公式，当p点位于物面s外部和远方控制面c的内部之空间域时，有如下公式：()由远方边界条件可知，远方封闭控制面上的积分趋于零，从而上式化为：()又由式〔〕可得：()得到混合分布模型，为了得到单一分布模型表示的扰动势，在物体内部域中构造一个适宜的内部解。于上述物体外部的点P，函数在物体内部域中没有奇点，在物体内部域中对函数和用格林第三公式，得：()将〔〕与〔〕相减，得：()取下式定解条件中的：(在内部)()(在)()那么式()成为：()其中，()2.2分布源密度的求解式()中右端分布源的法向导数极限由两局部组成，一局部是p点附近小曲面ε的奉献，另一局部是物面其余局部的奉献。法向指向取向物体内部，小曲面ε的奉献为2πσ〔p〕,那么有如下关系式：()再结合物面条件()，得到()这就是分布源密度所适合的线性积分方程。把积分方程()转换成线性代数方程组，即用离散量代替连续变量。把物面分成小块，记()用平面四边形或三角形来近似代替小曲面。具体做法如下，取第小块的四个顶点坐标之算术平均值，得到中心点的坐标。计算对角线向量的向量积〔指向与曲面法线指向相符合〕，用表示该方向上的单位向量，形成以为法线且通过中心点的平面，再把四个顶点向该平面作投影，以四个投影点为顶点组成平面四边形，用代替原来的小曲面，称为单元。通常把小范围内的分布源密度作为常数，因此只要分割不太粗，可以认为在单元上为常数，记作，从而()因此物面上的积分可以用个平面四边形〔三角形〕上积分之和来近似，即()上式左端的未知量是连续型变量，而上式右端的未知量是个离散量。为了求解这个未知数，须要个方程。取积分方程()中的动点为个单元的中心点，称之为控制点，即控制物面条件使之成立的点。用近似式(2.2.5)代替积分方程(2.2.2)的左端，便可以写出的阶线性代数方程组：()当计算出影响系数后，即可解线性方程组得到分布源密度。2.3速度势与附加质量的求解根据速度势在控制点pi处的值，由公式：()()根据2.2得到的分布源密度，求解线性方程组〔〕可得速度势的值。物体的附加质量，表示物体沿方向运动引起的方向的附加质量，公式如下:()()根据所求得的速度势的值可计算处附加质量的值。三、数值模型及参数计算3.1数值模型要求解流场中物体外表的速度势分布，需要先将物体的外外表进行网格划分。经过网格划分以后，原来的物体连续外外表被离散为N×M个相对独立的小平面，这些小平面构成了求解该问题的数值模型。Hess-Smith的根本思想是将连续曲面的积别离散为小单元来简化计算，其计算思路核心在于解该方程组：，通过求解线性方程得到。对于不同的计算目的，只需要改变控制面条件，即改变来实现。得到后，进而由求及附加质量，其中：为求，令，那么求得。故而，同理可得。3.2参数计算3.2.1计算对于球面：对于椭球面：设，，那么易知，，，故。〔a=1即为球〕对于圆柱面：在圆柱顶面上，同理在其底面上有，侧面上那么有3.2.2计算对于上述几何模型中的四边形上均匀分布奇点的诱导速度公式计算，首先将四边形上的分布源密度取为1，四边形四个顶点逆时针方向排列，顶点pi的坐标为〔ξi，ηi，0〕，其中坐标系oξηζ的原点为四边形形心，平面oξη即为四边形平面，那么需计算如下三个积分式：各积分式的计算公式为：在本几何模型中，需要将世界坐标系转换到局部坐标系中，查阅相关文献可知，当给定物体上的3个不共线的点，即可建立局部坐标系。其中局部坐标系的原点为点，轴的正向为的方向，轴的正向为的方向，轴的正向为的方向。设轴的单位矢量分别为它们在世界坐标系中的方向余弦分别用表示。对轴，可按下式给出：，，，为世界坐标系3根轴的单位矢量。矩阵就是由世界坐标系到物坐标系的坐标变换矩阵。经过坐标转换可得每个网格单元的Sx、Sy、Sz，然后得到世界坐标系下的S，即由于编程实现上述过程非常繁复，现采用一种更为简便方法的求解当时，当时，，从而求得的系数矩阵。3.2.3计算由可得当时，当时，如图，设位于上，取一小圆ε，转换为极坐标积分：〔广义积分收敛〕故可近似，得到系数矩阵。3.2.4计算m11、m33令，那么同理求得。对于双椭球体，如下列图所示，只需将第二个椭球的单元标记为〔，…〕即可。方程组的形式不变，在计算m11或m33时分别对两个椭球面的单元进行积分即可分别得到m11a、m11b、m33a、m33b。四、几何模型对于连续的积分方程，通常的数值处理是转换成线性代数方程组，即用离散量代替连续变量。在几何建模上，就是将物面s分成N小块，记为采用平面四边形或三角形近似代替小曲面。4.1球面网格划分4.1.1网格划分使用球坐标对球体进行网格划分，取球体半径R=1，以球心O为原点建立坐标系如下列图所示取xoz平面上的圆，圆上任意一点与x轴的夹角范围为0-π，划分N份，每相邻两份夹角为δφ=π/N，在球体上以球心为圆心，划分N条纬线。纬线为在球体上的同心圆，取任意同心圆，圆半径为r，将圆划分M份，每相邻两份之间夹角为δθ=2π/M,N个同心圆都划分M份，形成M条经线。由纵横交错的线将球体划分形成网格。取球体任意网格交点A,设坐标为A(x，y，z),由球坐标可将A坐标表示为：在matlab中根据以上网格划分原理划分网格，得到的球体几何模型为：4.1.2节点坐标取球体面上任意一个小网格，在小网格ABCD上，A标记为A(i，j)，即纵向第i个点，横向第j个点，那么i=1,2,…..,N+1j=1,2,……,M+1B,C,D可由A(i，j)表示，如上图所示。即可得到个小网格单元，故有个控制点P。对于顶点处，单元为三角形：或是4.1.3球面网格小单元面积对于四边形单元，AB//CD,由对称性可知，ABCD为等腰梯形。故ABCD为平面。ABCD的面积为δS。图解如图：梯形面积公式为：其中梯形高近似为腰长计算，得如下公式：EF的长度计算公式如下：梯形面积可表示为：4.2椭球面网格划分4.2.1网格划分同球形网面划分类似，得到N×M个小网格单元。与球面不同的是，沿x轴方向延长。由matlab划分网格得到的椭球几何图形如下列图所示。4.2.2节点坐标长短轴比为d。椭球面上任意网格交点坐标为:那么椭球体外表每个单元的控制点坐标为：4.2.3网格单元面积在O-XYZ坐标系中，椭球体的剖面图如下列图所示，易知：椭球面可由球面沿x轴延伸得到，即椭球面上某一小网格面ABCD的AC边延长为l2，由球面小网格面可知，面积为原来的l2/l1。即：即得小网格面积。4.2.4面积公式验证由3.2.3节得到了每个网格单元面积，由此可以得到椭球体外表积的近似值，为了减少计算过程中的误差，查阅相关文献后得到椭球体外表积计算公式，然后验证两者误差。设三椭球面方程的方程为：。其面积为：这里的，，又，分别为第一种和第二种椭球积分。对于，，k=0，那么F=u，E=u。当d=5，N=40，M=32时，计算出的单元面积和，理论值。相对误差为，误差很小，说明公式的正确性与精确性。4.3圆柱面网格划分4.3.1网格划分使用柱坐标对圆柱体进行网格划分，取圆柱底面半径R=1，以底面圆心O为原点建立坐标系。取yoz平面上的圆，圆上任意一点与y轴的所称角度范围为0-2π，划分K份，每相邻两份夹角为δθ=2π/N，在圆柱上按等距将圆柱沿高度方向分为M段，每份高δL那么由纵横交错的线可将圆柱侧外表划分形成网格。从而圆柱侧面可分为M*K个面元。按圆柱上下底面的极径均分为N份，每份长δr那么同理可以将两个底面分为2N*K个面元，总计将圆柱分为(2N+M)*K个面元。在matlab中根据以上网格划分原理划分网格，得到的柱体几何模型为：4.3.2节点坐标取圆柱任意网格交点A,设坐标为A(x，y，z),由柱坐标可将A坐标表示为：那么椭球体外表每个单元的控制点坐标为：4.3.3圆柱面网格小单元面积对于圆柱侧外表，AB//CD,AC//BD,易知ABCD为矩形。其中矩形沿圆柱高度方向边长为a=δL，另一边边长那么为对于圆柱顶面和底面可以按照两三角形面积之差计算五、计算结果显示5.1圆球附加质量计算结果经过不同组〔改变网格数量〕数值模型的求解，得到圆球附加质量系数M11、M33的计算结果，如下表4.1所示。表4.1单球体附加质量系数随球面网格数量的变化情况网格数50200450648968M110.60610.54550.52670.52110.5163M330.56290.52540.51620.51330.5108为了更加直观地考察附加质量随球面网格数量的变化情况，绘制出与上表4.1相应的折线图，这里将M11和M33的变化曲线画在同一张图里，便于观察比照，如下列图4.1所示。图4.1单个球体附加质量系数随网格数变化情况由上图4.1显示的曲线变化趋势可以看出，随着网格数的增加，数值计算的结果M11和M33都收敛于理论的解析解0.5。另外，所有的数值计算结果都大于理论值0.5。本项的计算结果，验证了球体外表在流场中的各向同性的性质，即对于球体，有M11=M33，且在一定程度上验证了解析解为0.5的正确性。或者反过来讲，数值计算的结果收敛于解析解，反映出HESS-SMITH方法的可行性和正确性。下面，从另一个角度来考察附加质量的性质。即改变物体模型，考察附加质量在各个方向上表现出来的性质如何。5.2椭球附加质量计算结果经过不同组〔改变网格数量〕数值模型的求解，得到椭球附加质量系数M11、M33的计算结果，如下表4.2所示。表4.2椭球附加质量系数随网格数量的变化情况网格数320490640980M110.06310.06410.06310.0608M330.92990.91540.93010.8942绘制出与上表4.2相应的折线图。这里因为M11和M33在量级上相差较大，故将M11和M33的变化曲线画在两张不同的图里，以便于观察各自的走向和变化。M11和M33的变化曲线分别如下列图4.2、图4.3所示。图4.2椭球附加质量系数M11随网格数变化情况图4.3椭球附加质量系数M33随网格数变化情况5.30双椭球附加质量计算结果对于计算两个椭球的干扰效应，计算时每个椭球划分为980个单元，轴间距按照短轴的3、5、7倍分别计算,得到两个椭球的附加质量系数。计算结果如表4.3所示：表4.3双椭球椭球附加质量系数随网格数量的变化情况间距比357M11a0.07610.06770.0644M11b0.07610.06680.0629M33a0.94890.86710.8493M33b0.94590.86300.8470依据表4.3绘制出图4.4及4.5图4.4双椭球M11随轴间距变化关系图4.5双椭球M33随轴间距变化关系5.4圆柱附加质量计算结果经过不同组〔改变网格数量〕数值模型的求解，得到圆球附加质量系数M11、M33的计算结果，如下表4.4所示。表4.4圆柱附加质量系数随球面网格数量的变化情况网格数375525750900M110.13420.13040.13460.1352M330.89820.90250.90500.9148绘制出与上表4.4相应的折线图。M11和M33的变化曲线分别如下列图4.6、图4.7所示。图4.6圆柱附加质量系数M11随网格数变化情况图4.7圆柱附加质量系数M33随网格数变化情况六、结果讨论误差分析6.1结果讨论对于单个球体，在无界流中其附加质量系数理论值是0.5，从计算结果及图表能看出，当球面单元划分数增加时，M11、M33均收敛于0.5。但M11、M33均大于0.5，这说明即便单元划分数目较大，但其外表仍不是光滑曲面，相比于球面，由水动力的知识易于理解其附加质量系数有微幅增量。图表知，当球面面元数为968时，M11、M33相对误差分别约为3.26%、2.16%，误差在允许范围内，初步验证了自编程序的正确性可可靠性。对于单个椭球，可从兰伯附加质量系数图表查得，当d=5时，其M11、M33理论值分别约为0.053、0.88,取面元数为980时的结果进行比拟，相对误差分别为14.72%、1.61%。M33的误差已经很小，可以认为计算结果相当准确。对于M11相对误差较大的原因稍后论述。对于双椭球，从图表能直观的看出两个椭球之间的扰动效应。当椭球附近流域内有其他物体时，其附加质量系数会增大，增量随间距为递减。理论上当两个椭球间距无穷远时，单个椭球的附加质量即为无界流中的附加质量，但从图中可以看出当轴间距为7倍以上时，单个椭球的附加质量就非常接近无扰动的值，此时可初步认为7倍间距为临界距离。另外理论分析知，两个椭球的附加质量系数应相等，从图4.4、图4.5可以直观了解到两个椭球的附加质量系数误差很小，与理论结果相符。对于圆柱，其M11的值未知，M33的理论值为1。从计算结果可以看出，随着网格划分数目的增加，M33的值逐渐趋近于1，当网格划分至900个面元时M33的相对误差为8.52%，与真实值比拟接近。推测其误差较大是因为网格划分过程中网格大小不均匀所致，位于圆柱上下底面圆心附近处的面元过于细长影响了运算结果的准确性。6.2误差原因分析1.方法误差〔截断误差〕在模型建立中，以小平面代曲面，本身就带来误差。求δS时，GH,EF的均舍去了高阶项O(x2),带来的极小误差。求Cij时，对于i=j,Cij≈2〔πδS〕1/2，由于是近似解，对于φ〔P〕的精准性有一定影响。此外模型建立过程中，网格划分如果长宽比过大，会使Cij的误差更为严重，而对于椭球而言，由于其结构较为细长，因此会产生大量细长的三角形单元；对于圆柱而言，在其底面上也会因网格划分过密导致这一问题，使得方法误差的影响加大。2.舍入误差在计算过程中由于计算机存储位数固定，数值在计算中会不停地在规定位数上取舍，带来舍入误差。即便计算时数值取双精度，但由于参与计算的参数数量大，计算过程长，造成误差积累，而求解过程中由于涉及到线性方程组的求解，因而舍入误差的量不容无视。由于M11量级相对较小，此时舍入误差占总误差的比值增大，因而此时会导致M11相对误差到达12%。6.3减小误差措施改变网格结构，使单元尽量接近正方形；适当增加网格数，增加方法近似性，减小截断误差。但数目不宜过大，会导致程序占用内存过大，计算速度明显下降，舍入误差累和增大等不利影响。参考文献[1]戴遗山,段文洋.船舶在波浪中运动的势流理论.国防工业出版社,2006.[2]范尚雍.船舶操纵性.国防工业出版社.1988[3]张亮，李云波.流体力学.哈尔滨工程大学出版社.2006.附录：matlab自编程序[1]圆〔椭〕球附加质量计算[2]圆柱附加质量计算[3]双椭球附加质量干扰计算[1]圆〔椭〕球附加质量计算N=input('输入N：');M=input('输入M：');R=1;%单元划分数N*M,球半径为1；d=input('输入d：');%长短轴比%划分结构化网格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;%单元划分数(2N+M)*K,圆柱底面半径为1；%划分结构化网格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;%单元划分数N*M,球半径为1；d=5;%长短轴比e=input('输入e：');%轴间距与短轴之比%划分结构化网格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+