周期性边界条件的一种处理方法

《周期性边界条件的一种处理方法》由会员分享，可在线阅读，更多相关《周期性边界条件的一种处理方法（10页珍藏版）》请在金锄头文库上搜索。

1、1953年3月数 值计算与计算机 应用第1期周期性边界条件的一种处理方法*黄艾香(西 安交 通大学)AME T H O DFO RT R E A T IN GPE R IOD IC A LBO UN DA R YCOND IT IO N SHuangAi一xi是犯ng(X乞赫J落a o忿。儿gU”公”rs公 才 y)AbstraCtTlli,paPergiv esan lethodforremovingperiod ie aleon strainteonditions undorthee onditionof maintainingfe aturesoft hesystemoffiniteel

2、ezn entalgebr aieequ atio n s.Thefeatureoftllismethod15thatthenun lberofvariable sandorde rofone一dimensioalar rangementremain sunehanged,sym metryandpo sitivedefiniteProPe rtiesar emain-七ainedandtheeholesk ymethodeanst illbeaPPlied.It15Pro vedthat th场metliod15extl,emelyeonv enientandeffeetinPraetie

3、e.一、引言在数学物理方法中,经常会遇到有部分边界满 足周期性条件的边值问题.例 如,在透平机械内部三元流动的流场延伸部分人工分界线(面)上(如图1中几2,几1 ),流函数沙满足 周期性边界条件叫r32训r, :+.又如 电机中的磁场分布(如图2),在 边界几2,几;半径相同处,其磁位A满 足周期性条件Ar3 :一一Ar,.若对具有这 类条件的微分方程运用有限元方法,那么在求解过程中得到 的有限元方程组就 必 须处理 周期性约束条件.在约束条件 是线 性情形时,一般 采用拉格 朗日乘子 法.这 种方 法需要增加求解参数.本文 给出一种在 不增加变量个数并保持方程组对 称正定 的情况下,解 除周

4、期性约束条件的 方 法.设有限 元方程组为Ax一B,其中(1)*1 98L年夕月2 6日收到.数 值计 算与 计算机应用1 9 8 3图1透平机械流面示意图、 、 、 . . . 古J了 户才甄甄l al aaa口月i口月2.召” ”图2电机转子截面示意图/bl 1bZI ”一(:b,户了才 口. . . . 砚 、一一AA为对称正定矩阵.在(l )式中,若某两个变量x ,勺有周期性约束x了=ax;一Q,其中a,Q 为常数.不失一般性,可设ij,我们 的目的是对方程组(l )解除约束(2)得方程组(2),使l )不改变方程的个数和变量的个数;2)保持 系数矩阵 的对称性与正 定性;3 )若(l

5、 )中系数矩阵是带状的,则尽量不 改变或少改变半带宽;斗)若(l )中系数矩阵是按一维排列顺序存贮,则不改变一维排列 顺序.二、解除周期性约束和矩阵变换设初等矩阵其中。i是第个分量 为1刃圣石,i(a)一其余分量 为0=l,刀,=I+(,一ei)e歹,(3夕的单位向量,(,二,纷)T.这里a,刀,=o(z铸i,了)对(l )作变换石,(a)Ax=石、j(a)B. (4)E,s(a)A一!E,i(a)才一才l知(4)与(l )是 等价的.用(2)代替(4)中第j个方程有Aix=Bi,(5)1期黄艾香:周期性边界条件的一种处理方 法其中、 、. . e e 1 | | | 布121 aijaliA

6、1一ao11a ii+aaiz二a ii+aa ia一,+aaaniB,一(b,b,一:,一口,bs+:,bi+abi,口月 月,b。)r.记其中向量妙,(井,终,杏,l,对(约作变换,杏:杏:二,I+(夸l,一ei)。歹,T.这里(6)、,产、/ . “t ) ”一(l+a11)a一a 、s,杏l,-一a ,i(l铸i,j)L、j(a)A:x一Lj(a)B,得AZx(7)其中AZL,(a)注;一(BZ,a汗硫 理碑 BZLi(a)B,这里姗一歼 奢一a卿一必彭a子,l,a子,(左铸i,=al,一 一a,j,友l,2,”,n)ac+a左i一aa左 saii+aaii,其余b*+a* s口(友粉

7、i,j,友一Q,占:2,=c夕+占i(l十az,)a十a1.a卿=a,=l,2,n),+abz,(8)必“( z )“卿) a.b.bJ, J . .e s l由于AZ匡,(a)刁:=1乙,(a)汉,一汉:,故(7)与(5)等价,而且是(4)解 除 了约束(2)后所得到 的方程组,也 是将(l )解除了 约束(2)的方程 组.从(s )可以看出,(7)仍然保持了系数矩阵的 对称性.下 面证 明A:仍为 正定矩阵.事实上,对A作合 同变换,变换矩阵取E i,(“),那么数 值 计 算 与 计算机应用1983自注(。=E,(a)AEz、(a)、龟龟. . |2 1 1 1 ! |a51+aas s

8、二ais+aassaii+aaii十aa is+aZass价,+“aa。了也是对称正定的.ani+aa”sa1 1aia la i+a na. /夕了l w e e e w e s e l l e e 、由于A(0 )与A是合同的,故A(0 )主 子式均大于零.由对称正定 的充要条件知,A(。的各阶对A:作合同变换,变换矩阵仍取凡i (a )得刁(2)一Es(a)AZEs,(a)a11+aaii口l”aii+aa iz+aaii+aZaz护,ai。十aasna。i+aaos:a- .aa(L r。.味a/犷了了| | ll|e e l . .、注意到A(2)与才(0 )仅第j行和第j列 的元素

9、不 同,主子式,得A(2)是对称正定矩阵,而A,与A(2)合同,A2)的各阶顺序主子式都是A(。的AZ也是对称正定的.而故三、解除周期性约束在计算机中实现的方法综上所述,为了将方程组(l )解除约束(2),并保持方程组的对称 正 定性,仅需按(8 )改变(l )中A的第j、i行和第j、i列 以及右端列 阵的 元素.由于有限元(线性)方程组的系数矩阵是对称正定,而且是带状的,因此在求解过程中,一般是按一维排列存贮下三角半带宽内元素.解除了约束(2)的方程组(7 )与(l )不但第i、j行有变化,它们的半带宽也有变动.因此我们必 须讨论 由(l ) 变换到(7)之后,一维排列改变的情况,以及其总长

10、度变化的计算方法.应注意(7)与(l )的其余各行的半带宽不变.设数组元素KD仅+l )为标定(l )中A第互行主 对角元素在一维排 列中的地址,并令K D(l )一。,乙*、二*分别表示第及行的半带宽和第一个非零元素的列 号,那 么有L*K D伏+l)一KD仅),二*=友一L*+l,左l,2,n.(9)再用I GP表示A中元素召*,(,*(l(劝在一维排 列中的地址,显然有IGp一K D(左+1)一寿+1.(10)l期黄艾香:周期性边界条件的一种处理方法由于(7)中第j行半带宽为1,故从方程组(l )变换到方程组(7 )时,第i行的改变皿 为15=Ls一1.(11)如果。,。,则第i行半带宽

11、不变,否则第行半带宽改变,其改变量为l、一,、一。s(i一z)一(乙一乙,).(12)因此,新系数矩阵的第萝行半带宽减少L i一1.第i行半带宽若增加的话,其改变量为l.显然,若l(乌一1,则第i行半带宽增加量 可 以被第i行半带宽减少量所补偿.否则系数矩阵的一维排列长度就要增加.由(11)和(12)知,其增加量 为S一l,一(Li一l)i一i一L+1.(3)第i、i行半带宽的改变,使得数组KD改变,矩阵A的一维排列长度NsKn(N+一)也要改变,其改变量由S所确定.实现方程组(l )解除周期性约束(2)而得方程组(力,其框图如下:l )变量及数组说明.NE周期性约束条件的个数,输人参数.N方

12、程组变量个数,输人参数.Ns一维排列 数组的长度,输人参数.DS系数矩阵一维排列长度的增加量,结果参数.I P(2,NE)存放(2)中每对自变量下标的数组,xP(l,NE)j,IP(2,NE)ti(ii),输人参数.AQ(2,NE)存放(2)中常数a、Q,AQ(l,NE)a,AQ(2,NE)Q,输人参数.A(Ns )方程组系数矩阵变半带宽一维排列数组,解除约束后为新方程组系数矩阵变半带宽一维排列数组.B(N)方程组右端项.KD(N十l )标定系数矩阵对角元素在数组A中的地址数组.一维排列增加量计算框图:SIC(NE,N+l,DS,IP,KD)D S0NEjIP(l,l)Ip(2,l)L,KD(

13、j+1)一Kn(j), 乙K o(i+1)一KD()数值计算与计 算机 应用19 83年l,nsj一Ls+l,*一乙i+115S=Ina X】11a X0,DS! o,一。i l z一(Li一l) l D S+S NE l R ETU RN解除周期性约束后A,B,KD的生成框图:PER(NE,Ns,Ds,N+l,N,IP,AQ,A,KD,B) S0 l自K=.0= 二一A E=c一N二二,一工=了乙,一二、A) ,+=IP(l,l),iIp(2,l)L,一KD(j+l)一K El),口AQ(2,l)。s一i一Ls+l,。;一乙l=maxo,。、mi 一左JGP=KD(j+l)一j+气IG p一

14、K D(i+l)一i+左B(左)+Q*x,A(JG P)一。 lV才A(JG p),B( 友)=ye s ,?一一飞 A x(左一。sl)左x*AIA(IGp)A(IGp)+x关A l 毒+|l1期黄艾香:周期性边界条件的一种处理方法产 . . . . . . .一-. . . . 口 . . .曰 . . .! 左一j一lyesAl(j一。,+1)4 一A 1,一j“i ?一几IGPKo(i+r)一i+j,+汉(IGp)左(K D(i+l)刀(KD(i+l)+(c+刀(IGP)*AI才(IGP)一AI召(i)刀(i)+B(i)*A I+c*l Q,B(I)=一Q,左(KD(z+1)一 l通!

15、+!+i一lm*=左一K D(左l)+K D(夜)+1用刃ye sKGpKD(左A(K GP)=+l)一友+z,x=才(K G p)o,B(交)B(友)+Q*xyes 及,?nO寺杏十 Al(友一m,l)一AI*x,IGpKD(i+l)一i一友才(IGp)理(IGP)+AI*x。*左一KD(左一)+KD(交)+l1m口yes!+lIG PKD伏十l)一友一i,JGp=KD认+l)一左一jx一才(JGp),B(友)B(友)+Q*x汉(xGp)=A(IGp)+AI*x,A(JGp)o 今 一一数值 计算与计 算机应用1 983年R ETURN四、应用实例根据上述方法,研制了标准子程序S IC、PER,它们分别用6912机和FO RTR ANw两种语言写成.从1977年开始,我们在计算各种 透平机械内部任意流面上流 函数的有限元解中,均使用这两个程序处理周期性 边界条件,取得了良好 的效果。透平机械内部三元流动任意流面上流函数价满足边值问题期黄艾香:周期性边界条件的一种处理方 法4 ,、l/d(必!;、:o,价: :一己ld 必bpo刀必】。3 :其中乙为流量,p试,句为密度.一沙r 3:+d,例 如,将此应 用于计算径流式透平压气机,流函数小在r3,、r32(图3)上满