《工程声学中的有限元法_沈崤(1)》由会员分享，可在线阅读，更多相关《工程声学中的有限元法_沈崤(1)（11页珍藏版）》请在金锄头文库上搜索。
1、声学第4期AC 丁A学报A CU ST IC 人195一年7月工 程声 学中的有限元 法*沈嚎孙洪生(中国科学院声学研究所)近 年来,有限元法已经在工程声学中应用.例 如,分析不规则房间内的声场,分析电声换能 器的振动问题以及用来评价抗性消声器的特性等.本文讨论用有限元法求解波动方程.用有限元法计算了不规则形状房间 的简正频率。比较了线性插值函数和三次插值函数的计算结果和本征值的收敛 性.通过二维模型的计算,证明了房间形状的不规则性对室内扩散声场没有影响.一、引言在声学现代化 过程中,广泛地应用 电子计算机,快速傅氏变换(FFT)是一个方面,有限 元法是另一个重要方面。有限元法是一种离散化 的
2、数值计算方法,它在理论推导中应用矩阵方法,实际计算则用 电子计算机.有限元法最初是在五十年代作为处理固体力学方法出现的,它是结构分析矩阵法的一个分支.初期的有限元法建立在虚功原理和最小位能原理上分析复杂结构的应力,19 60年以后分别用变分法和能量平衡法导出了有限元法解析方程组,使振动动态分析获得很大的进展.六十年代 中部分研究工作进人工程声学领域叩,例如分析不规则房间的简正振动方式【3 - -习,计算电声换能器的振动 6 1和消声器的传声特性即,等.本文讨论把波动方程转化为变分问题和用有限元法求解的基本方法,通过计算二维刚性封闭空间的简正频率的实例,具体说明插人函数的影响,不同方法收敛于准确
3、解的速度,最后给出不规则形状对小混响室简正频率影响的计算结果,证实与文献【 9 获得的结果一致.有限元法的适应性较强,是解复杂形状动态响应的有力工具,目前还需要进一步发展各种新的有限元模型,探索其数学理论基础,编制通用的计算程序,研究它在非线性声学中的应用.二、波动方程的变分间题声学中常常需要求解如下波动方程:vZP其边界条件有下列三种:户上旦 越c Za才2(1)fl(r)fZ(r)(2)+入户!,f3(r),力L nla一a力L 月,a一a乙气*1980年3月1 8日收到。声学学报19 8 1年式中p是 声压,是声速,为 边界,n为边界向外的法线方向,友、h是常数,fl,九,f3是边界上
4、的函数分布.式(2)是非齐次边界条件,它可以 通过适 当变换而化为齐次边界条件,因此 以后均讨论齐次边界条件,即式(2)的等号右边为零.对 于简谐振动p一P(,)e一沁,式(l )可写 为:v,户+左,P一0(3)式中功为角频率,友为波数,交。/ c.式(3)是亥姆霍兹(He l mh ohz )方程,用有限 元法求解这类问题时,首先把它变换为相应的变分问题.它等价于下列“商”泛函的变分问题:J (刃一D(刃/E(刃一临界值已砂(4)其中D(户)=(5)产1 1 .,丙. . , . . , .E(户)=;含,grd。,d V+,含。叨!;合”dF式中V是 边界包 围的区域.第二、第三类边值条
5、件自动为变分问题所满足,无须作为定解条件列出,称“自然边值条件”,第一类 边值条件应作为定解条件列出,称“强加边值条件气因此式(2)的三种边值条件仅剩下一种,即p,=0(6)这里已经变换成齐次边值条件.式(4一6)就是波动方程(l )的变分问题.它也可以直接从能量 平衡原理 导出,以 刚性边界 为例,封 闭空 间的动能和位能分别是:1Zp功21】gradP1UFPUV产恤 妞、,砂产 . 1万式中p为介质密度,为声速,。为角频率.在刚性边 界的封闭区域 中,自 由振动的能量守恒,即K一T一常数可 以看出上式是 与式(4)一式(5)等价的,因为是 刚性 边界,在 式(叼和 式(5)中h,0.三、
6、基本方法有限元法的基本方法是:把求解 区域分成有限个单元,在每个单元上构造近似函数,这个近似函数用结点坐标和结点参数来表示.因此,整个区域的变分问题就被离散成各单元 结点参数的求极值或临界值的问题,由于各个单元结 点是相互联系的,所以各单元的集合就成了整个区域 上对所有结点的求极值或临界值的问题,它 是一个齐次线性方程问题,用线性代数方法4期沈嚎等:工程声学中的有限元法求出问题的离散化 的近似解.归纳起来是以下六个步骤:(l )区域划分,(2)单元近似函数的构造,(3)单元分析,(4)总体合成,(5)边界条件的 处理,(6)最后用线性代数方法计算出结果.区域划分是按物理条件,把求解区域分成有限
7、个单元,单元的 形状、数量(大小)和分布是十分灵活的.但是,它们必 须遵守下述规则:单元的顶点相互连结,没有一个单元的 顶点是在其它单元的 边界上.实践表明还应使处 在边界上的单元的边 界尽量与物理边 界相吻合,使获得结果更精确、更符合物理 问题 的要求,整个区域中各个单元 的形状尽量选用 同一类型 以便于分析和计算.一般说来,单元数目越多,所得结果越准确.因此,在保证准确度的前提下,尽量减少单元的数量(即单元的尺寸大 些)以节省计算工作量.如果事先从物理分析出发能推断出区域中某些地方待求函数值变化乘缓或急剧 些,那末在平缓处可适 当减少单元数目,反之在急剧处增加单元数目,或在认为关键之处 增
8、加单元数目,以提高结果的准确度和相对地节省计算工作量通常在一维问题中采用线单元,二维问题中采用三角形单元、矩形单元,三维问题中采用四面体单元、长方体单元 等,如推广到多维问题那就是抽象的多维单元.为 了与物理条件吻合,有时也使用混合单元,例 如对于边被加厚的矩形膜片振动,其厚边用线单元,膜片用矩形单元.有时也使用 特殊形状的单元,例如扇形单元、柱面体单元等.但这样会给计算带来困难,因此通常不使用.显然,二维问 题中,三角形单元能够兼顾各方面的要求,实 际计算比较简单,因此最常用.本文讨论也 用三角形单元,如图1所示.。 邸. 一图1区域划分(习 和单元坐 标系统(b)D ividedr egi
9、onandeoo rdinatesysternfo relements单元近似函数通常采用多项式函数,这是因 为对它进行微积分运算比较方便,并且只要增加它的项数就 能以任意准确度逼近准确解.多项式的项数与单元 结点参数的个数相等.显然,若所选单元结 点参数越多,单元自由度越高,所求得的解 也越准确.对于二维问题来说,常用的 近似函数如下:线性插值函数为户一a,+。Zx+a3,(7) 二次插值函数为pal+aZx+a3y+a ;xZ+asxy+a6yZ(s)三次插值函数为声学学报1981年p一a :+aZx+a3y+a,xZ+asxy+a 6yZ+a 7x3+as(xZy+xyZ)+a gyZ(
10、9)它们分别有三、六、九个单元结点参数.单元结点参数通常选用单元结点上的函数值或 函数值的一次偏导数或 同时选 用这两种.单元结点则选在单元的顶点上,单元 边界上的特定位置(在二维问题中,单元 边界 上的特定位置 通常是边界线中点,在三维问题中则是边界面的几何中心),或者同时采用单元的顶点和单元边界上的特定位置.由于线性插值近似函数的运算最简便,虽然它的准确度较差,但可以用增加单元的数量来改进,所以最常使用它.计算表明,采用线性插值近似函数,收敛的速度较慢,因此二次插值近似函数和更高次插值近似函数的应用 也逐渐增加.实际上,以上二个步骤是把物理问题转化为数学问题.以二维刚性封闭区域简正频率的计
11、算为例,它相应的变分问题是式(4)和式(,),其中五.0.因为第二类边值条件被自动满足,不需列 出,所以边界条件的 处理就省去了.对于这样的 问题,用三角形单元来划分,并采用线性插值近似函数.因此,对单元e(如图1所示)来说:P(x,夕)a;。+a。x+a;,夕二,夕。(10)假定区域(图l )被分成N:个单元,共有N。个结点,对于每个单元有三个结点,结点即三角形的三个顶点.结点参数也 是三个,选用在三个结点坐标上的声压.把近似函数用结点坐标和结点参数来表示,并用矩阵形式写 出.令左()召发e)。;e)召;。T(一l)u(e)一“沪“犷)“;)T(12)z=lx夕(一 3)式中“沪,“沪,“沪
12、分别是 函数(声压)在单元结点l,2,3上的数值.字母右上方(。)表示单元。的值.式(10 )可 以写成:P(x,夕)二21A份,x,夕。(14)把结点坐标城。,诊,( i一1,2,3)及相应的函数值“:。代人上式可得:u 。万的才 心l(25)其中 “沪引 L H、”,一“ 飞 了 1二今,y飞JJ(1 6)仅与坐标有关.由式(z,)和式(14)可得.户(二,夕)z月(。一。 。x,夕。(17)上式即用结点坐标及其参数所表达的近似函数.Ni,(x,)一艺,份,(x,)对于整个区域来说,函数的近似表达式是:甲,(x,夕)x,夕G(18)斗期沈嚎等:工程声学中的有限元法式中甲。(x,夕)是“屋顶
13、”函数,甲(x,夕)lx,夕。0二,y 乡c而表示整个区域.把式(15 )代人式(4)和式(5),变分问题就离散成:J(户)J(“1,。2,帆)=D(夕)D(“,峋,。、。)(19)E(户)E(“,峋,u二。)式中D(一!。含!(豁)+(哥)“:一!。专尸“!(20)式(19)必须满足临界方程:O 二,.J又“1,“2,Ou idE,OD 泞.一U d“i“口“i,二二。)l,2,N 0(21)在式(1 9)一式(2 1)中,“*是整个区域上的结点参数,i1,2,N。,即把各单元的结点参数在整个区域内依次编号,1,2,N 0,它不同于单元结点参数,故去掉其右上角的(。).式(19)是变分问题式
14、(劝和式(劝的离散化方程,式(2 0)是逐个单元对近似函数进行具体运算,式(21 )是总体合成.在自然边值条件下,边值问题自然满足,无须考虑.最后便是对式(21)进行计算,可以证明它是一个线性齐次方程组.因此下面着重叙述式(2 0)的运算和式(21 )是怎样合成的。四、单元分析和总体合成单元分析就是计算式(2 0),为了便于最后总体合成,单元的三个结点的顺序是按逆时钟方向排列.在三角形单元,线性插值近似函数情况时,往往是一开始就把整个区域上的结点依次逐个编号,例如单元的三个结点编号是i,j,左(如图1所示).把式(17)代人式(20 )可以得到:刀。一生“(。T(犬l。+尺)r“。2石甸一生。
15、(。一T尺犷)z“(。12(22)声学1981年式中瓦气l 今乡“ b美J付纳举r. . . . . . .e s.L一一Kc声c毛c友;bib护对b天b了c ic矛c罗c ic友c夕c走(23)c*c,c又til矛t矛t走t萝z,t,又_r.l.I. . .L一一K2, t i转r. f.l. .L一一K其中b/一、义(y、一y搜,一众(,一,f,一戎(,一,小一众(“一,一众(一,、一众( x了一t ) ,(24)(25).介、一一片卜为八.l,l,钾一一气一1 2粼一6一一才考x矛x弋,土t工jl1一2 一一它即是单元e的面积.总体合成就是所有单元的集合,把式(2 2)代人 式(1 9),按矩阵叠加法则相加得:、少产O, 白了丫. 1, D(“,E(u l,u Z,u、。)一上“r了(尺11+r又21)“2“2,u、o)式 中一工“T尺3“I2一竺珠兰一丝2图2矩形闭室的例子I“一u,u Z,“、。丁Kl,KZ,【K3是N。阶方矩阵,它们分别是式(23)中尺l。,尺三亡),尺;)叠加而成,具体方法是:K矩阵中第i行,j列的元素是由相应的单元簇阵工K(勺 中,( e1,2,N口脚码为 l j 的元素的总和,脚码是指式(23)中元素的二个下标.如果所有单元中均 无脚码为 l i 的元素,则它为零.故二mpz。f:。eoangzarenezosodr。m今用 图2所示矩形闭室
