《复变函数论 Functions of a Complex Variable》课程教学资源（参考文献）二维弹塑性间题的复变函数边界元法.pdf_大学文库

第19卷第2期同济大学学报 Vol.1962 1991年6月 JOURNAL OF TONGJI UNIVERSITY Jun.1991 二维弹塑性问题的复变函数边界元法唐寿高 (工程力学系) 摘要本文将复变函数边界元发展用于解二维弹塑性力学问题，给出基于初应力荷载增量法的增量形式的弹塑性问题复变数边界积分方程及所有基本关系式，编制了相应的计算机程序并给出算例。与通常边界元解此类问题相比，本方法具有公式统一化，程序简单且多功能，计算机计算时间短且精度高等优点。计算表明，复变函数边界元特别适合用于需迭代计算的弹塑性等非线性问题，关键词：弹塑性，复变函数，边界积分方程。引言自从Swedlow和Cruseu1等人于1971年将边界元发展用于解弹塑性力学问题以来，许多学者对这一领域作了富有成效的研究，如文献[2]所述。就作者看来，其中具有代表性之一的是Banerjee等人的工作。他于l980年首次提出了基于初应力法的直接边界元法用于二维弹塑性分析31，所用初应力迭代步骤与Zienkiewicz]在有限元中所用方法相类似。后来 Banerjee和Raveendra2!引进加速初应力迭代技术和二次单元模式改进了原有的成果。最近，Henry和Banerjeets1通过数学处理得到的特别积分，避免了由于考虑非线性影响而要对区域内部计算的面积分或体积分。然而，上述方法仍有待于完善。如由于奇异积分，使解在边界及附近塑性区内难以收敛；考虑非线性因素后程序编制繁复，且无通用性，必须对所有边界划分单元，这样若采用迭代法求解，需耗费较多的计算机时间等等。本文的尝试为解决上述问题提供了一个可供选择的途径。文中仍采用惯用的初应力法，但由于核函数用复位势基本解来表示，所得增量形式的边界积分方程具有一般意义，所编程序可处理任意多个问题的基本解。针对具体问题采用相应的复位势基本解，这样可省去许多本来需要划分的边界单元，使方程组阶数大为降低，有利于反复迭代计算。再则解在不需分割单元的边界及近旁弹塑性区域内可给出令人满意的结果。还有本方法公式简单，程序编制容易且具统一性多功能，均比通常的边界元法优越。需要指出的是早在1982年Mukherjeel1就已将复变数用于边界元间接法分析蠕变和断裂问题。他也利用Muskhelishvilil)弹性力学复变数理论，但由于是直接从非齐次双调和方程变换得到边界积分方程，因而方程的意义及核函数等与文献[8]中的不同。此外，Mukher jee 的方法在实际运算时各主要量都是事先转变成实数来进行的。然而利用复变数的主要目的之本文收到日期：1990年3月6日 ?1994-2016 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net

第 1 9 卷第 2 期 1 9 9 1 年 6 月同济大学学报 J O U R N A L O F T O N O Jl t ) N I V E R S I T Y V o l . 19 恤 2 J U n . 19 9 1 二维弹塑性间题的复变函数边界元法唐寿高 ( 工程力学系 ) 摘要本文将复变函数边界元发展用于解二维弹塑性力学问题 , 给出基于初应力荷载增量法的增量形式的弹塑性问题复变数边界积分方程及所有基本关系式 , 编制了相应的计算机程序并给出算例。与通常边界元解此类问题相比 , 本方法具有公式统一化 , 程序简单且多功能 , 计算机计算时间短且精度高等优点。计算表明 . 复变函数边界元特别适合用于需迭代计算的弹塑性等非线性问题。关键词 : 弹塑性 , 复变函数 , 边界积分方程。一、引, . 言. -刁自从 S w e dl o w 和 C r u s e[ ` ’ 等人于 1 9 7 1 年将边界元发展用于解弹塑性力学问题以来 , 许多学者对这一领域作了富有成效的研究 , 如文献〔幻所述。就作者看来 , 其中具有代表性之一的是 B an e r je 等人的工作。他于 1 9 8。年首次提出了基于初应力法的直接边界元法用于二维弹塑性分析 s1[ , 所用初应力迭代步骤与 Z i e n ik e iw c z I 4〕在有限元中所用方法相类似。后来 B a en rj e e 和 R a ve e n d r a[ “ 1 引进加速初应力迭代技术和二次单元模式改进了原有的成果。最近 , H e n r y 和 B a n e rj e e[ ” 1通过数学处理得到的特别积分 , 避免了由于考虑非线性影响而要对区域内部计算的面积分或体积分。然而 , 上述方法仍有待于完善。如由于奇异积分 , 使解在边界及附近塑性区内难以收敛 ; 考虑非线性因素后程序编制繁复 , 且无通用性 , 必须对所有边界划分单元 , 这样若采用迭代法求解 , 需耗费较多的计算机时间等等。本文的尝试为解决上述问题提供了一个可供选择的途径。文中仍采用惯用的初应力法 , 但由于核函数用复位势基本解来表示 , 所得增量形式的边界积分方程具有一般意义 , 所编程序可处理任意多个问题的基本解。针对具体问题采用相应的复位势基本解 , 这样可省去许多本来需要划分的边界单元 , 使方程组阶数大为降低 , 有利于反复迭代计算。再则解在不需分割单元的边界及近旁弹塑性区域内可给出令人满意的结果。还有本方法公式简单 , 程序编制容易且具统一性多功能 , 均比通常的边界元法优越。需要指出的是早在19 8 2年M u k h er j e 「“ J就已将复变数用于边界元间接法分析蠕变和断裂问题。他也利用 M us k h e il s h vi il g l 弹性力学复变数理论 , 但由于是直接从非齐次双调和方程变换得到边界积分方程 , 因而方程的意义及核函数等与文献〔8] 中的不同。此外 , M uk h e r je 的方法在实际运算时各主要量都是事先转变成实数来进行的。然而利用复变数的主要目的之本文收到日期 : 10 90 年 3 月 6 日

188 同济大学学报第19卷一却是相同的，即在高度应力集中的边界上免去单元划分以便作精确的数值分析。Hromadka 等门将复变数用于边界元目的是对柯西积分进行离散，只解决了直接由柯西积分给出解的较为简单的标量场问题，而弹性力学等较为复杂的向量场问题的解是不能直接由柯西积分公式给出的。在边界元中使用复变函数基本解的概念首先是由黄弘量]提出，并对二维弹性问题作了成功的尝试。本文表明，复变函数边界元不仅适用于线弹性分析，而且特别适合于需迭代计算的弹塑性等非线性分析。二、复变函数边界积分方程及其基本公式基于加权残值理论的弹塑性力学平面问题增量形式的复变函数边界积分方程为 C(p)i(p)=」[U(p,qp.(q)-Pi(pq)i(q)]dr(g） +JnBp,g）m(gd0(g) (1) 式中 Cu(p)=-Pu(pq)dr(q) (2) 其中(g),(q)和°(q)分别是面力增量，位移增量和初应力增量，这里均作为权函数， U(,g》,P(中，q)和Bw(p,g)则是核函数，分别表示由于单位集中力在点卫沿x,方向作用，而在点9产生的沿：方向的位移、面力和应变（见图1）， x=1+x,i=√-I,以后当用作上下标时，i,i,k,1=1,2: 还有，初应力增量σ°(q)可由弹塑性总应变增量(q)通过本构关系表示为 (q)=(Di-D)m(q) (3) 式中，Di,为弹性本构张量，D?::为弹塑性本构张量，1(q) 为包含弹性和塑性应变的总应变增量。图1区城口及边界T示意通过边界积分方程(1)并考虑边界条件后可直接求出边界上 Fig.1 Diagram showing 所有单元上的(q)和(q)。于是，域内任意点的位移增量和应 region p with 力增量即可由(1)式计算得到，但此时C:(》=1,即 boundary r (p)=[U:(pq)p(q)-P(pq)i(q)]dr(q) +SBw(p(Qd(R (4) 考虑上式的微分，通过熟知的应力应变增量之间的关系可求得应力增量的积分表达式为 dp)=Dui(pa)p4(qdr(g）-JnSi(pq)in(gdr(g） +JTw(pq)on(qdQ(q (5) 其中应力影响系数 D4(p,q)=10U(p,q)+4[U,(p,g)+Uh,(p,q)] (6) Si(p,q）=161P1,(p,q)+4[P:,（p,q）+P1,(p,q)] (7) ?1994-2016 China Academic Journal Electronic Publishing House.All rights reserved. http://www.cnki.net

同济大学学报 1 9 卷一却是相同的 , 即在高度应力集中的边界上免去单元划分以便作精确的数值分析。 H r o m a d k a 等, ’ l将复变数用于边界元目的是对柯西积分进行离散 , 只解决了直接由柯西积分给出解的较为简单的标 t 场问题 , 而弹性力学等较为复杂的向量场问题的解是不能直接由柯西积分公式给出的。在边界元中使用复变函数基本解的概念首先是由黄弘量` “ l提出 , 并对二维弹性问题作了成功的尝试。本文表明 , 复变函数边界元不仅适用于线弹性分析 , 而且特别适合于需迭代计算的弹塑性等非线性分析。二、复变函数边界积分方程及其基本公式基于加权残值理论的弹塑性力学平面问题增量形式的复变函数边界积分方程为、尹. 1 1口0 吸矛夕、`、 e 一 ( , )“ ( , ) 二 J 厂〔。 ` , ( , , 。 ) , . ( ; ) 一 , ` , (。 , 。 )` · ( 。 )〕d r ( 。 ) + L 。 : 一 ( , , 。卜“ ` 。 ( 。 ) d 。 ( ; ) 式中 e 。 , ( , ) 二一 1 : , ` , ` , , 。 ) ` r ( 。 ) 其中矛 . ( 妙 , 公. ( 砂和云, ` 。 ( 妙分别是面力增量 , 位移增量和初应力增量 , 这里均作为权函数 , U ` , (乡一 ) , p 一 ( P , 口)和 B 一 , ( P , 口) 则是核函数 , 分别表示由于单位集中力在点 P 沿 x ` 方向作用 , 而在点 q 产生的沿幻方向的位移、面力和应变 ( 见图 1 ) , : = 劣: + x : i , i = 亿二丁 , 以后当用作上下标时 , i , j , k , l = z , 2 ; 还有 , 初应力增量云 . , 。 ( 妙可由弹塑性总应变增量云. , ( 妇通过本构关系表示为 d ` , o ( . ) = (刀了, 。一刀亨亨。 : ) · 云一: ( 。 ) ( 3 ) 式中 : D : ` 二为弹性本构张量 , D 丁亨 . , 为弹塑性本构张量 , 自, ( 妙为包含弹性和塑性应变的总应变增量。通过边界积分方程 (1 )并考虑边界条件后可直接求出边界上所有单元上的乡 ` ( 妙和 ` . ( 妙。于是 , 域内任意点的位移增量和应力增且即可由 ( 1) 式计算得到 , 但此时 C 。 , ( )P = 1 , 即图 1 F lg . 区域 O 及边界r 示意哪a召拍 m s五。 , 如. r铭i o n 0 w i th 加u n d a yt r ` , ( , ) = 丁 : 〔。 , , ( , , 。) , · (; ) 一 , · , ( , , ; 卜“ ( 。 ,〕` r (。 ) + 上 , ` . , ( , , ; ) · “ 一( ; ) d。 ( 。 ) ( 4 ) 考虑上式的微分 , 通过熟知的应力应变增量之间的关系可求得应力增量的积分表达式为 * ` , ( , ) 二【 _ D 二 , ( , , ; )矛。 ( 。 ) d: (。 ) 一份。。 , ( , , 。 )` . ( ; ) d r ( 。 ) J J , l , + J 。 : 一( , , , ,`一 ( 。 , d。 (。 ) ( 5 ) 其中应力影响系数 D一 J ( p , g ) = 浇6 . JU 一 , : ( P , g ) + 召[ U 二 , s ( P , q ) + U 一 , ` ( P , q )〕 S 一` , (乡 , q ) = 孟6 , J P 一 , , ( p , q ) + 拌[ P `一 ( p , g ) + P , 。 , . ( p , 红) 〕 ( 6 ) ( 7 )

第2期唐寿高：二维弹塑性问题的复变函数边界元法 189 T1i(p,9)=16:1B:,(p,9)+[B:i,(p,q）+B1Ai:(p,q)] (8) E 这里()：=0)/2w0z/ax。式中1=1+-2，k=21+70对于平面应力问照，只须将弹性模量石换成+E,泊松比？换成！中，即可。三、用复位势基本解表示的核函数根据弹性力学复变函数理论「8，)，位移核函数可用复位势表示成 U)(3-4)tnRe((z +δ2lm(P,(zzp)门-[ò1Re(zp,(zzp) +:(22))-δ121m(2p,(zz)+:(xz)]} (9) 其中z=x1-x2i;p(zzp)和(zzp)为复位势基本解，()'=d()/dz。由于单位力在点，沿1方向作用而引起点x的应力可从下式求得 1:(z,z)=6:Re[2p1(zzp)+（-1)(zp1(zzp)+1(zz)] +(δ62+6,26i1)Im[zp(z,zp)+1(z,zp)] (10 于是，容易求得面力核函数 P1i（2,zp）=a01ia(2,2） =nd:Re[2p(z,zp)+（-1)(zp(zzp）+0:(xzp)门 +（n162+26:1)Im[zp1(zz）+1(zz)] (11) 其中为边界法线”关于轴x。的方向余弦。最后，应变核函数B1(z,z)用复位势表示如下： Bm22w）=2[U,2）+Uw()门 =3-4[3Rc(e2w）+p.22w》 +δ2Im(Pi(zzp)+p,(z2p))门 -[OmRe(z(i(z,z)+ii(zz))+i(zzp) +.(x,zp)）-6Im(z(P,（z,z)+9.（zz) +.i(z,2p）+1,(zzp))]} (12) 至此，二维弹塑性力学问题的复变函数边界元所有基本公式已全部给出。这些表达式即可用来编制相应的计算机程序。由于这些公式并未涉及复位势基本解的具体内容，因此可用来分析具有各种特定边界的弹塑性平面问题，只要导出相应的复位势基本解并编制计算这些基本解及其导数的子程序以供调用。与通常边界元只能处理一种基本解相比，功能大为加强，这是因为通常边界元需用基本解的具体表达式来编制程序，如采用不同问题的基本解则整个程序须重新编制。再则这里所有表达式最终由一对复位势P和中来表达，这就使得基本公式最大程度地简化，程序的编制与调试也就很方便了。 ?1994-2016 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net

第 2期唐寿高 : 二维弹塑性问题的复变函数边界元法 T 。: ` s ( P , 口) = 又6 ` , B : 一: , l ( P , 口) + 产〔B : * 。 , , ( P , 口) + 刀 ; 。 , , 、 ( P , 口)〕 ( 8 ) 这里 ( ) , ` = a( ) / 0勺。加可ax ` 。式中沐 = E梦 ( 1 + , ) ( 1 一 2 , ) ’ 升 = E 2 ( 1 + v ) 。对于平面应力问题 , 只须将弹性模量 E 换成 ~ 具牛黑洲、 1 十 7 , - E , 泊松比 , 换成1 + F 即可。三、用复位势基本解表示的核函数根据弹性力学复变函数理论 [ 8 , 。 ] , 位移核函数可用复位势表示成 U · , `: , : , , = 命 “ 卜 `“ )〔` j I R。 ` 甲` : , : · , , + 6 , : I m (甲` ( z , z , ) )〕一〔占, , R 。 (牙甲飞( : , : , ) + 叻` ( : , z , ) ) 一己, : I m (三甲奋( : , z , ) + 叻` ( : , z , ) )〕} 其中牙二 x : 一劣: i ; 甲( : , : , ) 和叻( : , z , )为复位势基本解 , ( ) ` = d ( ) / d z 。由于单位力在点 z , 沿 l 方向作用而引起点 z 的应力可从下式求得口: ` , ( z , : 一 ) = 饥, R 。〔2甲飞( z , z , ) + ( 一 i ) ` (牙中飞( : , : , ) + 叻飞( : , 言 , ) )〕 + ( 6 : : 6 , : + 6 . 2 6 ` , ) I二〔牙中叮( z , : , ) + 功飞( z , z , ) 〕于是 , 容易求得面力核函数 P I ` ( 之 , 2 , ) = 作。口 t i 。 ( z , 之少) 二 , j 6 ` , R 。〔2中; ( : , z , ) + ( 一 1 ) ` ( 牙中叮( : , : , ) + 叻: ( z , z 一 ) )〕 + ( 称: 6 。: + ” : 6 ` : ) I 优〔三甲丫( z , z , ) + 劝几( z , 之一 ) 〕其中 : . 为边界法线 ” 关于轴 x 。的方向余弦。最后 , 应变核函 , 数 B , , 。 ( : , 介 )用复位势表示如下 : ( 9 ) ( 1 0 ) ( 1 1 ) B ` J。 (之 , z 户) = 1 , , 丫丁 L U ` 一八 z ’ 之一 ) + U ,一 ( 之 , 2 , )〕 1 , , ~ 二一厂 ~ 弓气石一 4 科 4拼) [ 6 一: R 。 ( 甲` , , ( : , z , ) + 甲 , , ` ( : , : , ) ) + 己。 Z l m ( 甲` . , ( z , z , ) + 中 , , ` ( : , : , ) )〕一〔6 一, R 。 (三(甲 ; , , ( z , : , ) + 甲 ,i , ` ( z , z , ) ) + 叻` , , ( : , z , ) + 叻, . ` ( z , : , ) ) 一占. : I m (三(甲飞 , , ( : , ’z ) + 甲 ,I . ` ( z , z 一 ) ) + 协` , , ( : , : , ) + 叻, , , ( : , : , ) )〕} ( 1 2 ) 至此 , 二维弹塑性力学问题的复变函数边界元所有基本公式已全部给出。这些表达式即可用来编制相应的计算机程序。由于这些公式并未涉及复位势基本解的具体内容 , 因此可用来分析具有各种特定边界的弹塑性平面问题 , 只要导出相应的复位势基本解并编制计算这些基本解及其导数的子程序以供调用。与通常边界元只能处理一种基本解相比 ,功能大为加强 , 这是因为通常边界元需用基本解的具体表达式来编制程序 , 如采用不同问题的基本解则整个程序须重新编制。再则这里所有表达式最终由一对复位势甲和叻来表达 , 这就使得基本公式最大程度地简化 , 程序的编制与调试也就很方便了

190 同济大学学报节19卷四、复位势基本解这里只列出两种问题，即无限平面和开圆孔的无限平面问题的复位势基本解。然而读者可以文献[8,9]所述方法导出其它边界问题相应的基本解，例如椭圆孔，方孔及半无限平面等形式的孔口或直线边界问题。计算时不需在这些边界上划分单元（如局部边界有约束或力作用，只需在这约束段进行单元分割)，这是因为基本解事先已满足了相应的边界条件。 1.无限平面基本解（图2) 考虑了位移对称性后，此时的复位势基本解为町 p1(z3)=A1n(z-z） (zzp)=B1ln(2-zo)-A1z/z-z）+aj》 (13) 式中A,=1/儿8(1-4)]；B1=(-1)'(3-4)·A;i=√/-1， 1=1,2,是荷载作用点，2是欲计算函数值的点（下同）。 2。开有圆孔的无限平面基本解图2无限平面示意 Fig.2 Notation for 此解p9和可以通过在无限平面基本解(13)上叠加一对关 infinite plane 于圆孔的映像函数p:和：得到「8)，即 P9(zzp)=PI(zzp)+pi(z3zp) (14) 功9(zzp)=1(220)+1(232p) 其中p,(z,zp)和1(zzp)即为(13)式所示，映像函数p(z,zp),:(z,z®)在考虑了位移对称性后成为 pi(zzp)=B1n[z/八z-z)门-A1(zp-za)za/[7(z-z)门 (zzp)=A血[z/八z-2a门+A(Z。-R/z)/(z-za) (15) -B:R[1/z-1/(z-zo)3/z-A(zp-z)z3/[z(z-z.)2]-A 其中R为圆孔半径，。=R/为p关于圆周的像点。我们已看到，在无限平面的解上叠加了使之满足圆孔自由边界条件的映像函数即得到问题2的解。其它边界问题的复位势基本解也类似处理8)。五、方法的实施在程序编制中采用了有限元边界元常用的基于初应力法图3开圆孔的无限平面示意的增量荷载方法。将待求问题的有关边界进行单元划分，并 Fig.3 Notation for infinite 对估计会屈服的区域布上适当大小和数量的三角形网格，考 plane with hole 虑了边界条件后，增量形式的边界积分方程(1)成为 [A们{8}={i}+[B]G°} (16) ?1994-2016 China Academic Journal Electronic Publishing House.All rights reserved. http://www.cnki.net

同济大学学于设第 19 卷四、复位势基本解这里只列出两种问题 , 即无限平面和开圆孔的无限平面问题的复位势基本解。然而读者可以文献 8[ , 刘所述方法导出其它边界问题相应的基本解 , 例如椭圆孔 , 方孔及半无限平面等形式的孔口或直线边界问题。计算时不需在这些边界上划分单元 ( 如局部边界有约束或力作用 , 只需在这约束段进行单元分割 ) , 这是因为基本解事先已满足了相应的边界条件。 1 . 无限平面若本解 ( 圈 2 ) 考虑了位移对称性后 , 此时的复位势基本解为 , “ 二 ] 么一胜甲一 ( 之 , z 一 ) = A I ln ( 之一 z , ) 护. ( z , : , ) = B , ln ( 之一 : , ) 一 A ,牙, / ( z 一 : , ) + 诬,} ` 1 3 , 式中 A : 二 i “ + 1 , /〔s 兀 ( z 一科) ] ,刀: = ( 一 i ) ` ( 3 一 4拼) · A : ; i = 了二丁 , , = 1 ,2 , z , 是荷载作用点 , : 是欲计算函数值的点 ( 下同 ) 。 2 . 开有国孔的无限平面签本解、夕. úì J `了 Jl. 、、卜l . J 此解叫和叻“ 可以通过在无限平面基本解 ( 1 3) 上叠加一对关于圆孔的映像函数叫和必得到『’ , 叼 , 即中? ( z , z , ) 二甲 , ( z , z , ) + 甲丁( z , z , ) 叻, ( z , z , ) = 功: ( z , z 一) + 功丁( z , 之一 ) 其中甲. ( 之 , z , )和叻: (之 , z , ) 即为 ( 1 3 )式所示 , 映像函数叫 ( 之 ,介 ) , 性后成为中 , ( z , z , ) = 石: 玩〔 z / ( z 一 z 。 )〕一亘: ( z , 一 z 。 ) z a / [牙, ( : 一 z 。 ) ] 吵1( z , z , ) = 迈 : ln [ : / ( z 一 z 。 ) 〕+ A : ( 玄。一 R Z / 之 ) / ( : 一之。 ) 图 2 无限平面示意 F i s . 2 N o t a ti o n 企o r i n f i n i t e lP a n 。 ( 1 4 ) 川 ( 之 , : , ) 在考虑了位移对称一石: R 至 [ 1 / : 一 1 / ( : 一 z 。 )〕/ z 一迈 : ( : , 一 z 。 ) 之盆/ [ : ( z 一 z J ) “ 〕一 A : 其中 R 为圆孔半径 , 之。 = R V几为 z , 关于圆周的像点。我们已看到 , 在无限平面的解上叠加了使之满足圆孔自由边界条件的映像函数即得到问题 2 的解。其它边界问题的复位势基本解也类似处理『8 , ” ] 。五、方法的实施图。开圆孔的无限平面示意 F 宜: 3 No 扭t i o n fo r I n f i n i t e 川a n e w i t卜 h of e 在程序编制中采用了有限元边界元常用的基于初应力法的增量荷载方法。将待求问题的有关边界进行单元划分 , 并对估计会屈服的区域布上适当大小和数量的三角形网格 , 考虑了边界条件后 , 增量形式的边界积分方程 ( 1) 成为 [ A 〕{充} = {乙卜 [丑〕笼行 o } ( 1 6 )

第2期唐寿高：二维弹塑性问题的复变函数边界元法 19t 式中[A]为2N×2N非对称满阵，}为2N×1的待求量列阵，其元素为边界上的未知位移增量和未知面力增量，{}是边界上的已知值与矩阵[G]和[F]中相应列的乘积的2N×1列 ◆ 阵，[B]为2N×3M矩阵，这里M是已屈服的三角形网格数目，{⊙}为3M×1列阵，其元素为已屈服的三角形网格内的初应力。计算时首先在全部外荷载作用下，作为弹性问题初算，得出应力{a}和位移{4，并令 L=amax/c,这里G为等效应力，0，为材料的屈服应力。我们乘用Mises.屈服准则。若L ≤1则属于弹性问题，计算结束。否则当前应力{σ，}和位移{u.}可分别通过如下二式求出并贮存 {ge}={ap}/L{u.}={u}/L 然后将超弹性极限部分的荷载P(工一1)分成适当的荷载步长以通常边界元所述方法迭代计算。但本文用如下方法处理奇异积分：以该单元（或域内三角形网格）的中点作为新的结点，将单元（或网格）分成二个子单元（或三个子三角形网格），再分别进行数值计算。结果表明，这样的处理方法是方便和有效的。本文所用数值积分公式为Gaussian积分公式和Gauss- Radan积分公式I,)。若用一次或高次元模式对程序进行修改，还可进一步提高计算精度。六、计算实例及讨论这里就弹性和弹塑性问题给出两个算例。例1设一无限平面开有一半径为的圆孔，孔内对径受压（图4a)。取集中力P= 1,000kg/cm,r=3cm,作为平面应变问题。我们只在荷载作用处各布置了一个单元（图4b)。除这二处以外，其余孔边及近旁均可求得与解析解很为吻合的计算结果（见表1）。 4元② 例2一开孔铝质板条见图5中所示。板条二端受均布拉力巢度为σ，材料硬化元系数H=2206.496MPa,屈服应力o=238.302 MPa,=0.2,E=68646.55MPa。作为标准例子，此问题曾被Theocaris和Marketos1l, (a) (6) 图4对径受压圆孔 3 Zienkiewicz,及Banerjee()等人分别作过 Fig.4 Infinite plane with hole subjected 实验及有限元和边界元等的详细分析。这里为 to a couple of force P 便于计算中间薄弱截面上的应力，我们考虑对称性，仅取上半板条分析而不取1/4板条。计算结果与所提其他方法结果的比较示于图5和图6中。从上述两例可看出本方法与通常边界元相比，由于复位势基本解在域内和圆孔边界上都是精确满足的，而通常边界元基本解只在域内精确满足，故无论是弹性区域或塑性区域，特别是孔边及其近旁，本方法均可得到更为精确的计算结果。第一例仅用二个单元，不到通常边界元所需单元的1/10，CPU时间约为通常边界元所需的1/8。第二例本方法无须在孔边划分单元，这不仅使得方程组的阶数大为降低，节省了大量的计算机时间，而且在该边界及近 ?1994-2016 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net

第会期唐寿高 : 二维弹塑性问题的复变函数边界元法式中〔A〕为 ZN x ZN 非对称满阵 , 硬戈}为 ZN x l 的待求量列阵 , 其元素为边界上的未知位移增量和未知面力增量 , {朴是边界上的已知值与矩阵〔 G 〕和〔 F 〕中相应列的乘积的 ZN x l 列阵 ; 〔 B 〕为 ZN x 3 M 矩阵 , 这里 M是已屈服的三角形网格数目 , {旬为 3M x l 列阵 , 其元素为已屈服的三角形网格内的初应力。计算时首先在全部外荷载作用下 , 作为弹性问题初算 , 得出应力 {` } 和位移{勒} , 并令 L 二于。。 /久 , 这里子为等效应力 , 久为材料的屈服应力。我们乘用 M i s e s 屈服准则。若 L ( 1 则属于弹性问题 , 计算结束。否则当前应力 { .a }和位移 { : 。 } 可分别通过如下二式求出并贮存 {吼 } = {价 } / L { 。。 } 二咬, , } / L 然后将超弹性极限部分的荷载 P ( L 一 1) 分成适当的荷载步长以通常边界元所述方法迭代计算。但本文用如下方法处理奇异积分 : 以该单元 ( 或域内三角形网格 ) 的中点作为新的结点 , 将单元 ( 或网格 ) 分成二个子单元 ( 或三个子三角形网格 ) , 再分别进行数值计算。结果表明 , 这样的处理方法是方便和有效的。本文所用数值积分公式为 G au s s ia n 积分公式和 G a us s - R a d a n 积分公式 I” ,川。若用一次或高次元模式对程序进行修改 , 还可进一步提高计算精度。六、计算实例及讨论这里就弹性和弹塑性问题给出两个算例。例 1 设一无限平面开有一半径为 7 的圆孔 , 孔内对径受压 ( 图 a4 ) 。取集中力 P = 1 , o 0 k g/ c m , , = c3 m , 作为平面应变问题。我们只在荷载作用处各布置了一个单元 (图 4石) 。除这二处以外 , 其余孔边及近旁均可求得与解析解很为吻合的计算结果 ( 见表 1 ) 。例 2 一开孔铝质板条见图 5 中所示。板条二端受均布拉力集度为口。 , 材料硬化系数 H = 2 2 0 6 . 4 9 6 M P a ,屈服应力乡。 = 2 3 8 . 3 0 2 M P a 声 = 0 . 2 , E = 6 s 64 6 . 5 5 M P a 。作为标准例子 , 此问题曾被 T h e o 。 a r i s 和 M a r k e t o s t” ] , Z i e n k i e w i c z [`了, 及 B a o e r j e e [ , J 等人分别作过实验及有限元和边界元等的详细分析。这里为图 4 对径受压圆孔 Pi g . 4 】n f i n i t e P l a n e w 皿t h h o l o s u b j e e t e d t o a e o u P l e o f fo r e e P 便于计算中间薄弱截面上的应力 , 我们考虑对称性 , 仅取上半板条分析而不取 1/ 4 板条。计算结果与所提其他方法结果的比较示于图 5 和图 6 中。从上述两例可看出本方法与通常边界元相比 , 由于复位势基本解在域内和圆孔边界上都是精确满足的 ,而通常边界元基本解只在域内精确满足 , 故无论是弹性区域或塑性区域 , 特别是孔边及其近旁 , 本方法均可得到更为精确的计算结果。第一例仅用二个单元 , 不到通常边界元所需单元的 1 / 10 , C P U 时间约为通常边界元所需的 1 / 8 。第二例本方法无须在孔边划分单元 , 这不仅使得方程组的阶数大为降低 , 节省了大量的计算机时间 , 而且在该边界及近

192 同济大学学报第19卷旁可得到令人满意的计算成果，克服了通常边界元呈现的发散现象。表1X1轴上孔边及附近正应力与解析解比较 ×3=0 本文方法解析解 =/m 01(98.088/Pa) 0(98.066/Pa) 01(98.066/Pm) 0x(98.066/Pa) 5 0.1×10- 212.268 0 212.207 3.25 15.024 180.288 15.019 180.238 3.5 25.677 154.123 25.668 154.083 3.75 33,096 132.569 33.085 132.540 4.0 38.115 114.708 88.102 114.687 4.5 43.264 87.348 43.250 87.338 6.0 44.422 67.938 44.409 67.935 8.0 30.954 21.367 30.953 21.372 Tab.1 Stress comparison of numeric and analytic solutions on and near the boundary of x:axis 关于用初应力法分析弹塑性问题时要对域内屈服区域的“残余应力”进行体积分，目前国内外已提出多种方法来消去这些区域积分。作者认为，在本方法中采用类似处理原则上是可行的。 0.5 i 0.4 1.a 3 g 1.0 60.3 实验值I0] 旦0.8 实验值[10] 0.6 Aa边界元[3] 0.2 一般边界元[3] 本文方法 --…有限元法[4] 0.4 。一本文方法 0.2L 1.0(01.2 1.41.6 1.82.0(') 0.5i.01.52.02.3.035 图6 Een/oo 薄弱戴面AA'上的正应力结果比校（当2% 图5第一屈服点A处最大应变结果比较 =0.91时) Fig.5 Comparison of maximum strain Fig.6 Comparison of stresses at section at first yield point A determined AA of perforated plate determined by different methods by different methods 七、小结发展了基于初应力荷载增量法用于弹塑性力学分析的复变函数边界元直接方法。本方法程序简单，具有可扩充性，计算精确可靠，适合分析有高度应力集中的开孔弹塑性问题。由于一般只需较少的边界单元，故复变函数边界元适合于分析需迭代计算的弹塑性等物理非线性问题，特别是各种开孔开口问题，只要导出相应的复位势基本解即可利用本文所编程序作弹塑性数值分析。本方法的局限性是只能分析二维问题，这是因为复变数难以表达三维问题所致。 ?1994-2016 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net

第 2 期唐寿高 : 二维弹塑性问题的复变函数边界元法 19 3 本文主要内容为作者于 1 9 87 一 1 9 8 9 年初 , 由中英友好奖学金的赞助在英国曼彻斯特理工学院进修时科研工作的一部分 . 得到导师 R . K i t c h i n g 教授的热情关怀 , 支持和有益的讨论 . 稿件整理时还得到同济大学创业基金部分资助 . 在此一并深表谢意。参考文献 [ 1 ] S w e d l o w J L , e t a l . F o r m u l a t i o n o f b o u n d a r y i n t e g r a l e q u a t i o n s f o r 3 D e l a s t o p ! a s t i e f l o w . I n t . J . S o l id s , 1 9 7 1 ; 7 : 1 6 7 3 一 1 68 3 . [ 2 〕 B a n e r j e e P K , R a v e e n d r a 5 T . A d v a n e e d B E M a n a l y s i s o f t w o a n d t h r e e d i m e n t i o n a l p r o b l e m s o f e l a , t o · P l a s t i e i t y . I n t . J . N u m . M e t h o d s E n g n g , 1 9 8 6 ; 2 3 : 9 8 5一1 0 0 2 . [ 3 ] B a n e r j e e P K , C a t h i e D N . A di r e e t f o r m u l a t i o n a n d n u m e r i e a l i m p l e · m e n t o f t h e BE M f o r t w o 一 d i m e n t i o n p r o b l e m s o f e l a s t o 一 Pl a s t i e i t y . I n t . J . M e e h . S e i . , 1 9 8 0 ; 2 2 : 2 3 3一 2 4 5 . [ 4 ] Zi e n k i e w i e z 0 C . T h e f i n i t e e l e m e n t m e t h o d , L o n d o n : M e G r a w 一 H i ll , 19 7 7 : 4 5 0一4 7 1 . [ 5 〕 H e n r y D P , e t a l . A n e w B E M f o r m u l a t i o n f o r t w o a n d t h r e e d i m e n - t i o n a l e l a s t o p l a s t i e i t y u s i n g p a r t i e u l a r i n t e g r a l s . I n t . J . N u m . M e t h o d s E n g n g , 1 9 8 8; 2 6 : 2 0 7 9一2 0 9 6 . [ 6 ] M u k h e r j e e 5 . B o u n d a r y e l e m e n t m e t h o d s i n e r e e p a n d f r a e t u r e , L o n d o n : A P Pl i e d S e i e n e e P u b li s h e r s , 1 9 8 2 . [ 7 〕 H r o m a d k a 1 T V , G u y m o n G L . A e o m p l e x v a r i a b l e b o u n d a r y e l e m e n t m e t h o d : D e v e l o p m e n t , I n t . J . N u m . M e t h o d s E n g n g , 1 9 8 4 ; 2 0 : 2 5一3 7 . 〔 8 〕黄弘量 , 唐寿高 . 复变函数论的边界元方法 . 同济大学学报 , 19 8 6 ; 1 4( 2) : 1 7 9一1 9 1 。 [ 9 〕 M u s k h e l i , h v i li N 1 . 数学弹性力学的几个基本问题 , 北京 : 科学出版社 , 1 9 5 8 - 〔1 0〕 S t r o u d A H . G a u s s i a n q u a d r a t u r e f o r m u l a s , C lif f s , N e w J e r s e y : P r e n t i c e - H a l l , I n e . E n g l e w o o d , 1 9 6 6 . [ 1 1 ] S t r o u d A H . A p p r o x i m a t e e a l e u l a t i o n o f m u l t i p l e i n t e g r a l s . C li f f s , N e w J e r s e y ; P r e n t i e e 一 H a ll , I n e . E n g l e w o o d , 1 9 7 1 . [ 1 2〕 T h e o e a r i s P 5 . M a r k e t o s E . E l a s t o 一 p l a s t i e a n a ly s i s o f p e r f o r a t e d t h i n s t r i p s o f a s t r a i n h a r d e n i n g m a t e r i a l . J . M e e h . P h y s . S o l id s , 1 9 6 4 ; 1 2 : 3 7 7 一 3 9 0 . 汗芦及