第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 。 第二例本方 法无须在 孔 边 划 分单 元 , 这不仅使得 方程组 的阶数大为 降低 , 节省 了大量 的计算机时 间 , 而 且在该边 界及近