内法线方向怎么求 曲线的内外法线向量怎么求

摘要：结构在静流场中振动时,附加质量的分布形式对其振动特性有非常大的影响,尤其是对较为轻薄的结构,但截至目前,对附加质量显式分布的研究甚少。基于三维问题的边界元法,在合适的Dirichlet及Neumann条件下,提出了一种较简单的方法,可求解结构以任意给定模式在不可压缩单相静流场中振动时的附加质量显式分布。数值计算表明,该方法的求解结果与理论解及实验结果吻合良好,可较好地应用于复杂形状结构在有界或无界静流场中振动时的附加质量分布计算。

关键词： 力学 数值计算 结构振动 附加质量 静流场

当一个物体在流体中作刚体运动或模态振动时,会带动其周围的流体一起发生运动,其中运动流体的动能可以用相应的附加质量进行度量。若假定流体无黏无旋,则依据势流理论,完全浸没物体的附加质量取决于其形状及运动模式,而与位移大小、速度及加速度无关[1]。

对于附加质量,除极个别简单情形存在解析解外,一般需由数值方法或实验进行计算。基于CFD(ComputationFluidDynamics)方法得到的力与加速度关系,Uchiyama[2]和Benaouicha等[3]分别求解了二相流中振动圆柱以及空泡流中翼型的附加质量,Fu等[4]利用离散涡方法研究楔形体入水时的附加质量矩阵,而Liu等[5]则通过对实验结果进行参数识别得到附加质量。

当较为轻薄的结构在流体中作模态振动时,附加质量会对其自振特性产生很大影响。基于动压力的显式表达,Han等[6]针对流体中的振动柱壳求解自振频率和振型,并提出了附加质量的简易表达式。Minami[7]采用涡片法求解作用于薄膜上的压力,并依据能量原理求得附加质量比。孙旭峰等[8]利用奇点配置法推导了三维情形下结构振动诱导流场附加质量的计算式。Yadykin等[9]则对一悬臂板以前10阶弦向模态振动时的附加质量进行了讨论。

上述研究所针对的均为总体附加质量,而实际上当结构在流体中发生模态振动时,附加质量的分布形式对其振动特性也有很大的影响。Li等[10]曾对一圆形膜片提出过附加质量的分布形式并进行了实验验证,但截至目前,附加质量分布的理论研究还极少。基于边界元法和一定的Dirichlet及Neumann条件,本文提出了一种较为简便的数值方法,可用于求解复杂边界条件下静流场中振动结构的附加质量分布。

1、振动结构附加质量分布的数值算法

1.1 三维Laplace方程的边界元解法

对于三维势流问题,控制方程可写为

式中,ϕ(x)为域V内关于任意点x的势函数。在域边界S上,Dirichlet条件和Neumann条件分别为

式中:S=Sϕ∪Sg;n为域V在边界S上的外法线方向。

由Green等式,奇异点P的解可写为

式中:q为边界S上的场点;ϕS(P;q)为式(1)的基本解

式中,r(P,q)为P和q的距离。

使P在边界S上趋于点p,即可得到边界积分方程:

式中,Sδ为包含点p且半径为δ的球面,故对光滑边界有c(p)=1/2,而对角点则有c(p)=φ/4π,其中φ为角点处域内一侧的内角弧度数。

将边界划分为Ne个单元,并以插值函数将边界变量ϕ和∂ϕ∂n∂ϕ∂n进行离散化,就有一组关于结点未知量的线性代数方程组。由于该方程组的系数矩阵一般为满阵,所以当自由度超过10000时,使用高斯消去法等直接解法将非常困难。此时,可采用快速多极边界元法[11,12],该算法可将运算数量级从N3降至N。

1.2 附加质量及其分布

设结构的前r阶模态为Ψ1,Ψ2,…,Ψr,则位移和速度的正则表达式可分别写为

式中,ξi为正则坐标。

若结构的振幅足够小,则其引起的周边流体运动可视为无旋流动,从而可写出势表达式

ϕ=∑i=1rϕidξidtϕ=∑i=1rϕidξidt(10)

式中,ϕi为对应于第i阶模态的单位势。

设结构边界SS被划分为NS个单元,则ϕi可写为

ϕi=∑j=1NSϕjiϕi=∑j=1ΝSϕij(11)

这里ϕjiij为第j单元的势,并满足

∂ϕji∂nj=−∂ϕji∂n=ψji∂ϕij∂nj=-∂ϕij∂n=ψij(12)

式中:nj为j单元指向流体域的内法线方向;n为外法线方向;而Ψjiij则为j单元的模态速度,其正负号由nj确定。

这样,式(10)中由第i模态分量引起的周边流体动能可写为

K=−ρ2∫SS(ϕidξidt)∂(ϕidξidt)∂ndS≈−ρ2(dξidt)2∑j=1NSϕji∂ϕji∂nAj=ρ2(dξidt)2∑j=1NSϕjiψjiAj (13)Κ=-ρ2∫SS(ϕidξidt)∂(ϕidξidt)∂ndS≈-ρ2(dξidt)2∑j=1ΝSϕij∂ϕij∂nAj=ρ2(dξidt)2∑j=1ΝSϕijψijAj (13)

式中:ρ为流体密度;Aj为j单元的面积。

注意到流体动能亦可写为

K=12mi(dξidt)2Κ=12mi(dξidt)2(14)

则由式(13)、式(14)可得

mi=∑j=1NSmji=ρ∑j=1NSϕjiψjiAjmi=∑j=1ΝSmij=ρ∑j=1ΝSϕijψijAj(15)

式中:mi为对应于第i阶模态的正则附加质量;mjiij为j单元的附加质量,其所代表的即为结构边界上的附加质量分布。

2、数值算例

2.1 以单位速度在无限流体域中运动的单位半径球

该问题附加质量的经典解为m1=23πρm1=23πρ,球面边界上的最大速度势则为ϕSmax=0.5。用边界元法求解该问题时,可用足够大的立方体模拟无穷远边界,无穷远边界条件取Dirichlet条件ϕ¯=0ϕ¯=0或Neumann条件g¯=0g¯=0,而球面边界条件取Neumann条件g¯=−vn=−cosθg¯=-vn=-cosθ,如图1所示。

随球面划分单元数NS的增加,得到的附加质量m1、球面最大速度势ϕSmax如表1所示。由此可知,随网格划分数量增加,计算收敛于精确解,而由式(15),NS=12744时的附加质量分布如图2所示(取ρ=1,下同)。

图1球面边界上的Neumann条件

图2单位半径球在流体中以单位速度平动时的

球面附加质量分布(NS=12744)

2.2 流体中平动和转动的无限长漂浮棱柱体

文献[13]给出了流体中的漂浮棱柱体以单位速度平动和转动时附加质量的半解析解,其中棱柱的转动中心为C,如图3所示。用边界元法求解时,边界条件如下:S1上ϕ¯=0,S2ϕ¯=0,S2上ϕ¯=0ϕ¯=0或g¯=0,S3g¯=0,S3上g¯=0g¯=0,其余边界条件如表2所示。设棱柱体的宽度为2(即a=1),长度取20,湿周划分为2922个单元,则其竖向平动、水平向平动及转动时的单位长度附加质量计算结果如表2所示(需要注意的是,表2中的附加质量结果已经消除了边界效应的影响,其余算例也作了相同处理)。

图3流体中平动和转动的无限长漂浮棱柱体

由表2可知,无论棱柱体下方的流域为无限还是有限,采用了合适边界条件的边界元法求解结果与理论解都非常一致。图4为e=∞时,采用三维边界元法求出的各种运动模式下的附加质量分布。

图4漂浮棱柱体的附加质量分布(e=∞)

2.3 沿径向振动的圆柱壳

设无限流域中的单位半径圆柱壳以径向速度cos2θ发生振动,如图5所示。则有边界条件

g¯=∂ϕ∂r∣∣r=1=−cos2θg¯=∂ϕ∂r|r=1=-cos2θ(16)

图5沿径向振动的圆柱壳(r=1)

设势函数为如下形式

ϕ(r,θ)=F(r)cos2θ(17)

并将其代入式(1),则有

∂2ϕ∂r2+1r∂ϕ∂r+1r2∂2ϕ∂θ2=0∂2ϕ∂r2+1r∂ϕ∂r+1r2∂2ϕ∂θ2=0(18)

这样就有

F⋅⋅(r)+1rF˙(r)−4r2F(r)=0F⋅⋅(r)+1rF˙(r)-4r2F(r)=0(19)

求解式(19),并由式(16)所示的边界条件就可以得到

F(r)=12r2F(r)=12r2(20)

根据势流理论,振动引起的周边流体动能为

K=−ρ2∮r=1ϕ(r,θ)∂ϕ(r,θ)∂rrdθ=14ρπΚ=-ρ2∮r=1ϕ(r,θ)∂ϕ(r,θ)∂rrdθ=14ρπ(21)

从而可知,该问题的附加质量理论解为0.5ρπ。采用三维边界元法求解该问题时,若取圆柱壳长度为20并将柱面划分为3206个单元,则可求得单位长度上的附加质量为1.56846ρ,误差为-0.15%,附加质量分布如图6所示。

图6圆柱壳表面的附加质量分布

3、静止空气中振动的膜结构

3.1 以基本模态振动的一维平面膜

图7所示为一两边固定的无限长平面膜结构。设该膜以半正弦基本模态振动,则在图示坐标系下,位移可表示为

h(x,t)=asin(πxl)cos(ωt),0≤x≤lh(x,t)=asin(πxl)cos(ωt),0≤x≤l(22)

式中:a为最大振幅;ω为圆频率;l为膜宽。

图7以基本模态振动的二维平面膜结构

基于薄翼理论,Minami曾推导出该膜结构的附加质量为0.68ρl,文献[14]则采用声学理论求解出附加质量为2πρl≈0.64ρl2πρl≈0.64ρl。采用三维边界元法求解时,膜的上表面和下表面边界条件分别为g¯=−sin(πxl)g¯=-sin(πxl)和g¯=sin(πxl)g¯=sin(πxl)。若取膜宽为l=4,膜长为20,并将膜的上、下表面分别划分为2076个单元,则可以得到附加质量为每单元面积1.74405ρ。

这个结果比Chen等的研究结果小了约31.5%,比Minami等的研究结果则要小35.9%。但是当我们重新审视Chen等的研究时,可以发现其推导结果实际上是单侧膜面与流体接触时的附加质量为1πρl1πρl,而双侧膜面与流体接触时只是简单地乘了2。Minami等在推导附加质量时,并没有直接求出膜片上、下的压力差,而是间接采用了涡片强度基本解γm。若在三维边界元法求解时对膜上、下侧的边界条件都取为g¯=−sin(πxl),g¯=-sin(πxl),则可得结果为2.5725ρ,与Chen等的研究误差为0.5%,与Minami等的研究误差则为-5.4%。由此可见,Minami等和Chen等的推导过高估计了膜结构的附加质量。

在上述尺寸和单元划分条件下,该一维平面膜的上表面或下表面附加质量分布,如图8所示。

3.2 静止空气中振动的圆形平面膜

Li等依据不同气压下自振频率的变化,对图9所示的周边固定圆形平面膜进行了附加质量测定,其中实验膜面的直径为300mm并可施加不同的预应力水平。由于测点布置的缺陷,导致该实验中部分三阶、四阶模态数据缺失,故本文仅以一阶、二阶模态的实测平均值作为比较。

图8膜的上(下)表面附加质量分布

图9Li等研究中的圆形平面膜实验装置

图10为圆形平面膜的前两阶振型,其中膜面划分为1866个单元。在以三维边界元法求解时,采用归一化模态速度作为膜的上表面和下表面的Neumann边界条件。

图10周边固定圆形平面膜的前两阶模态

边界元法的附加质量求解结果如表3所示。表中同时还列出了Li等的理论解以及一阶、二阶模态的实验平均值。需要说明的是,从Li等给出的位移能谱密度曲线可知,一阶模态的谱峰并不显著,而二阶模态的谱峰能量则要远高于一阶模态,故此应该说,二阶模态的附加质量实验识别结果可信度较高,实验误差相对较小。由表3可知,边界元法的求解结果与实验值非常一致,尤其是可信度较高的二阶模态附加质量。同时,在上述单元划分条件下,附加质量在膜的上表面(或下表面)的分布如图11所示。

4、结论

附加质量及其分布形式会在很大程度上影响静流场中振动结构的振动特性。但截至目前,对附加质量的研究主要还是集中在总体附加质量,而很少涉及其分布形式,对流场中发生模态振动的结构则更是如此。本文基于三维边界元法和经典势流理论,在适当的Dirichlet和Neumann边界条件下提出了一种可以求解静流场中附加质量及其分布的数值计算方法。通过典型算例的比较分析,可以得到如下结论:

(1)随着单元划分的细化,本文数值结果收敛于精确理论解,说明该方法是准确可靠的。

(2)不论是对无界静流场还是有界静流场,亦或是结构发生刚体运动还是模态振动,该方法得到的总体附加质量都与经典理论解吻合良好,而静止空气中振动膜结构的附加质量实测结果也与本文方法所求得的结果非常一致。

(3)虽然本文在推导过程中采用了模态速度及其正则坐标,但实际上,Neumann边界条件可以适用于静流场中以任意模式产生振动的结构,而不仅仅局限于刚体运动或模态振动。

就附加质量的计算而言,经典势流理论仅能给出一些简单情形的解析解,而有限元或有限体积等数值方法则只能以隐含形式间接考虑附加质量效应的影响,无法获得其显式分布。本文所提出的方法通过直接指定边界条件,可求得以任意模式在无界或有界静流场中振动结构的附加质量及其分布,对相关工程应用而言具有实用参考价值。