（固体力学专业论文）界面端奇异行为的数值分析.pdf

浙江大学博士学位论文( 2 0 0 2 ) 摘要 首先，从多重应力奇异性叠加的奇异应力场出发，通过将奇异性次数分离， 提出了一种利用常规的数值分析结果来确定多重应力奇异性( 包括振荡应力奇异 性) 次数的数值分析方法，并且对双材料以及三材料界面端模型的平面变形和三维 变形问题进行了数值计算。 其次，基于横观各向同性压电材料空间轴对称变形问题的通解，利用特征展 开法，给出了奇异点附近的位移场和奇异应力场。在此基础上，对具有任意界面 角和结合角的横观各向同性双压电材料空间轴对称界面端一般模型的轴对称变形 问题进行了理论分析，给出了该模型界面端的奇异性特征方程以及界面端附近的 位移场、电势、奇异应力场和奇异电位移场。为了检验理论分析的正确性，应用 有限元程序对轴对称界面端的奇异性次数进行了数值求解。 最后，从最小势能原理出发，在仅仅考虑奇异性支配区域这一前提下，对于 弹性接合材料的平面变形问题和拟平面应变问题，以奇异点为原点分别建立极坐 标系和圆柱坐标系，通过分部积分消除r 项，从而使奇异性问题的求解由原来的二 维降为一维：对于三维变形问题，以奇异点为原点建立球坐标系，通过分部积分 消除，项，从而使奇异性问题的求解由原来的三维降为二维。然后，采用有限元方 法，把降阶后的区域离散，得到了确定奇异性的特征方程。用数值方法求解特征 方程，就可以得到满足条件o r e n ) o 这时，4 和爿：都是实数，即 仁一学 ，爿，一! 二塑二丝 ( 2 l o ) 2 d 对( 25 ) 式取对数，得 1 0 9 0 盯l 一爿二盯：1 ) = l o g lb 1 1 一 l o g 产 i 0 9 0 盯一爿口：1 ) = l o d b ：l 一 1 0 9 ； ( 2 一l 】) 由上式可知，在应力奇异性的支配区域内，双对数分布l 0 9 0 q 爿：盯：i ) l o g i 和 1 0 9 q q 一爿盯：1 ) 1 0 9 i 的斜率即分别对应于两个不等的实应力奇异性次数 ，如。 3 、62 4 臼c o 这时，爿。和一：是一对共轭复数，即 爿，：互：一。+ 叫，：一宰 c z 一- z ， 由上式可知，在应力奇异性的支配区域内产生振荡应力奇异性。 ，毋，瓦，日。， 。和兀p ) 都是复数，代入( 2 5 ) 式得到 ：b p ”；一( 1 )( 2 1 3 ) 式中 盯：瓜石再百万，占：厢 y 。= t a n “【爿，。：b 一彳。a ：) 】，。= t a n 。( 毋) ( 2 1 4 ) 列( 21 3 ) 式取对数，可得 l o g 盯= l o g b 一 r l o g 声t 。l o g p = 8 l o g e 一丑，l o g 产 ( 2 1 5 ) 因此，若爿。，爿，己知，在奇异应力场的支配区域内，双对数分布l o g 口l o g ；和 。l o g e l o g 为两条直线，其斜率分别为振荡奇异性次数a 的实部厶和虚部2 ，。 1 2 浙江大学博士学位论文( 2 0 0 2 ) 在实际计算爿，和一，( 或对于复数情况为。和一，) 时，可在同一方向臼= 岛上， 以比例常数p 选取多个按等比级数分布的点，由连续的三个点组( t ，七+ 1 ，+ 2 ) 可 分别求出一系列爿。和爿，。在奇异应力场的支配区域内的同一方向上，爿和一! 应 为常数。因此，由最小二乘法可确定出爿。和爿：。 2 2 2 三重戍力奇异性 对于具有三个应力奇异性次数的异材界面端或界面角点，在以奇异点为原点 的极坐标系( r 臼) 中，奇异点附近的应力场可表示为 咔垡掣+ 茎掣+ 茎掣 ( 2 _ 1 6 ) 。 r o 。r 如r 式中- 、k ，和，p ) ( = 1 ，2 ，3 ) 分别为应力奇异性次数、与应力强度系数有关的系 数和角函数。引入无量纲因子i = r 三后，( 2 一1 6 ) 式可写成 e = 箸，( 臼) + 告厶( 臼) + 箸厶p ) ( 2 ) 为了确定- ，在某一给定的方向雠，) 上，任意选择三个应力分量巳，代入 到( 21 7 ) 中可得 声一 声七 ； h ，= h ，慨) ( i ，= l ，2 ，3 ) h = 2 l2 22 3 h 1 11 2h 1 j 【h ， 见2h 3 3 j ( 2 1 8 ) ( 2 一1 9 ) ( 2 2 0 ) 1 3 如如一巩 h = t，【，j q q 巳 ，j、，l 中式 浙江大学博士学位论文( 2 0 0 2 ) ，一 声一2 2 # 一 巨 其中，c ，= 日门h ，b 。是矩阵h 的元素片。的代数余子式。 在奇异应力场的支配区域内按等比级数选取四点，与奇异点的距离分别为 ，j = p “1 _ 心= 1 1 3 4 ) ( p 为比例常数) ，对应点的三个应力分量分别记为 一1 女= l ，2 1 3 ；j = l ，2 ，3 ，4 ) ( 其值可由普通的有限元法或边界元法等数值方法确定) ， 代入( 2 2 1 ) 式，得 等等织：t 瓦习i i 耳百刀刮， 矧善搿， 百孑而再瓦万一 矧等搿咄 ( 2 _ z z ) c ，盯f 4 ) + c 2 ，一4 ) + c ，盯叫o 式中 t ，= p o = l ，2 ，3 ) ( 2 2 3 ) 由( 2 2 2 ) 式可得 b ”一t ，口j 2 ) c ，+ b 一 ，仃1 2 1 ) c ：，+ b ! ”一t ，盯f 2 1 ) c ，= o b f ”一k ，盯 j ) c 。+ b p 一，盯p ) c ：，+ b ”一i ，盯1 3 ) c ，= o b f 3 一女，盯f 4 ) ) c ，+ b 5 ”一t ，盯! ) c ：，+ b p 一七，盯5 4 ) ) c ”= o ( 2 2 4 ) 要使关于c 。c ! ，c 、，的齐次线性方程组有非零解，系数矩阵行列式的值必须为零。 将行列式展开r 可得关于女的一元三次方程式 武，3 + 6 t ，2 + c 女，+ d = o ( 2 2 5 ) 式中 口= 盯f 4 b 1 3 ) 盯! ”一d 2 叫3 ) + a 2 毋一d 3 ) 硝2 ) + o 4 3 ) 盯! ”一盯f 2 ) 仃1 3 ) 1 4 、l，q q h 儿 岛 浙江大学博士学位论文( 2 0 0 2 ) 6 = d 4 b ! ”。p 一。c 3 ，。+ a ，b p d ”一0 1 ) 仃1 3 ) + 毋妍j 仃p l o ，以1 ) c = d 4 埘一础盯+ 砖1 d ”一盯f 2 q 1 ，) + 毋斜2 ) 盯 ) 一卅，盯笋，) d = c r f ”( 以盯l ”一吐2 d 1 ) 十叫3 2 q ”一一1 以2 1 ) + d 3 1 以”一a f ：q ) ) ( 2 2 6 ) 对方程式( 2 2 5 ) 求解时，存在两种可能的情况： l 、若( ，= l ，2 ，3 ) 为实数，则i 个实应力奇异性为 五，= l o g 七，l o g p ( 22 7 ) 2 、若t 。为实数，女：和屯为一对共轭复数，即实应力奇异性和振荡应力奇异性 并存。女：可表示为 足2 = t = 七r + 腑= p 厶+ 胁= p 厶【c o s ( 丑l n p ) + j s i n 0 ，l n p ) 】 ( 2 2 8 ) 即 置r = p 厶c o s ( 五，l n p ) ，七，= p 厶s i n n l n p ) ( 2 2 9 ) 于是，实应力奇异性为 五= l o g 七l l o g p( 2 3 0 ) 振荡应力奇异性的实部和虚部分别为 厶= l o g 如；+ 女? ) ( 2 l o g p ) = t a j l 1 0 ，k ) l n p ( 2 3 1 ) 2 2 3 多重应力奇异性 对于具有多重应力奇异性的异材界面端或界面角点 表示为 口：争生型 z 一 j = 1 r 。 参照上一节中的方法，可以得到个一元次方程 奇异点附近的应力场可 ( 23 2 ) 乜 七，“十d 一l 七，”+ + d l t 十d 。= o ( 2 3 3 ) 求解方程( 2 3 3 ) ，即可得到全部的奇异性次数。需要指出的是，采用同样的方法， 也可以利用位移分量来确定奇异性次数。 1 5 浙江大学博士学位论文( 2 0 0 2 ) 2 3 数值计算实例 为了检验上述分析方法的有效性，本节对平面应变状态下的双材料模型以及 三维异材界面端模型的应力奇异性次数进行了数值求解。平面模型的数值计算利 用通用的二维边界元程序b e m 2 d f l 6 0 】，而三维模型的数值计算采用通用的有限元 程序a n s y s 。 2 3 1两个实奇异性的二维界面端模型 双材料界面端模型如图2 1 所示，材料常数为：e = 1 0 8 g j p a ，v = o 3 3 e 二= 3 0 4 g 尸口，v ! = 0 2 7 。卅点附近应力奇异性次数的理论解为 ，= 0 4 9 1 4 7 9 五，= 0 2 l8 6 5 2 。数值计算利用通用的二维边界元程序b e m 2 d ，采用二次等参元。 界面端4 附近的单元网格划分如图2 2 所示。在单元划分时，预先按p = 2 布置有 关节点。爿点附近界面上的应力盯和r 。的双对数分布如图2 3 所示。由图2 _ 3 可 知，双对数分布为两条相交的曲线，这意味着在选用不同的应力分量、不同区间 的近似直线的斜率来确定应力奇异性次数时，将得到不同的结果。这表明该界面 端的应力奇异性次数不是一个，已不能直接用求解只有一个奇异性时的双对数应 力分布来确定界面端爿附近的应力奇异性次数。利用( 2 1 0 ) 式求得多组爿，和以， 其分布如图2 _ 4 所示，由最小二乘法可求得爿，= 1 3 9 8 3 4 3 ，以= 2 _ 3 8 5 5 7 9 。然后 根据( 2 一1 1 ) 式作出两条双对数分布1 0 9 华1 0 9 参，如图2 5 所示，两条 双对数分布具有很好的直线性，其斜率即分别为应力奇异性次数 ，= 0 4 9 0 5 15 ， 五，= 0 2 2 2 6 4 2 ，与理论解相比，相对误差分别为o 2 ，1 8 。较弱的应力奇异性 的误差相对较大，但仍有很高的精度。 浙江大学博士学位论文( 2 0 0 2 ) y 2 i ， 尹 岛”2 。 彳 2 图2 1 计算模型 墨 宝 j 飞 占 鼍 旦 l o g ( r ) 图2 3 a ，和t 。的双对数分布 ， l o 科叫) 图2 5 ( q 一爿h ) 的双对数分布 飞 i 图2 2 界面端剧驸返的网格划分 最小单元尺磋= 2 o 1 0 。7 缈 l o “，矿) 图2 4 系数爿。的分散性 f ， p l ! ! 。! 图2 6 计算模型 1 7 一q手吾一 t引牟叫十引士 浙江大学博士学位论文( 2 0 0 2 ) 2 3 2 振荡脊异性的二维晃面端模型 双材料界面端模型如图2 6 所示，材料常数为：e = 1 0 0 0 g 强，y 。= o 3 ； b = l g p a ，y 2 = o - 2a 以点附近的振荡奇异性次数的理论解为靠= o 3 0 5 8 5 9 ， 丑= o 1 0 3 2 0 3 。界面端爿附近的单元网格划分如图2 7 所示。在单元划分时，预先 按p = 2 布好有关节点。一点附近界面上的应力q 和的双对数分布如图2 8 所 示。由图2 ，8 可知，双对数分布明显为两条曲线，利用( 2 一1 2 ) 式求得的多组4 。和彳， 如图2 9 所示。由图中可见，以和爿，的分散性都很小，由最小二乘法可求得 如= 1 - 4 9 4 5 0 5 ，4 ，2o 8 7 7 3 5 3 。然后根据( 2 1 5 ) 式作出1 0 9 盯l o g 参和 l 。g p l 。g 专两条分布曲线，如图2 - 1 0 和图2 1 l 所示。两条分布曲线具有很好 的直线性，其斜率即分别为振荡奇异性次数a 的实部和虚部，即厶= o 3 0 6 5 1 5 ， 五，= 0 1 0 4 1 3 2 ，与理论解相比，相对误差分别为o 2 和o 9 。 图2 7 界面端爿附近的网格划分 最小单元尺寸o 2 1 0 。7 矿 i o 猷r ) 图2 8 q 和勺的双对数分布 - 1 8 浙江大学博士学位论文( 2 0 0 2 ) l o g ( ，) 图29 系数爿，和爿，的分散性 旦 墨 l o g ( r ) 图2 1 0 盯的双对数分布 l o 甙r ，h 图2 1 1 p 。的单对数分布 2 3 3两重实奇异性的三维异材界面端模型 双材料界面端模型如图2 1 2 所示，立方体l 的材料常数为：e ，= 3 9 0 g p a p ，= o 3 ；立方体2 的材料常数为e = 1 4 0 g p a ，v ，= 0 4 。数值计算利用通用有限 元程序a n s y s ，图中的模型被划分成1 3 5 4 1 个四面体单元，共有2 0 1 8 4 个节点， 最小的单元尺寸为5 l o ，计算时问大约为3 0 分钟。在单元划分时，由于无法 预先按一定的比例常数在立方体内的某一方向上布置有关节点，因此数值分析所 需的应力分量盯，和盯，只能通过插值得到。界面端d 点附近沿直线0 h 上的应力口。 和盯的双对数分布如图21 3 所示。由图21 3 可知，双对数分布为两条相交的曲线， 这意味着在选用不同的应力分量、不同区间的近似直线的斜率来确定应力奇异性 1 9 3 5 2 5 5 0 2 0 浙江大学博士学位论文( 2 0 0 2 ) 次数时，将得到不同的结果，即界面端的应力奇异性次数是多重的。由第四章中 的特殊有限元法可得到，界面端有两重实应力奇异性，奇异性次数分别为 丑= 0 2 5 5 8 ，t = 0 4 4 3 4 。将插值得到的应力分量代入( 2 一l o ) 式，可以求得多组爿。 和一，其分布如图2 1 4 所示，由最小二乘法可得4 = 0 0 1 3 1 6 ，爿，= 1 2 4 6 8 7 。然 后根据( 2 1 1 ) 式，作出两条_ 双对数分布j 0 9 0 巳一。盯： p ) 】o g 口) 如图2 1 5 所 示，两条双对数分布具有很好的直线性，其斜率即分别为应力奇异性次数 丑= o2 4 3 5 ，丑= o 4 6 】8 。与由特殊有限元法得到的数值结果相比较，相对误差分 别为48 和4 2 。 卜矿十百叫 告 罾 图21 2 三维界面端奇异性分析模型图213 盯。和盯：的双对数分布 l o g ( ，口) 图2 1 4 系数一的分散性 6 t 葡 旦 1 0 9 ( “d ) 图2 1 5 j ，一4 一。哎的双对数分布 2 0 - 浙江大学博士学位论文( 2 0 0 2 ) 2 3 4三重实奇异性的三维异材界面端模型 双材料界面端模型如图2 1 6 所示，材料常数为：e l = 1 0 4 g p a ，v = o 3 ： e 2 = 2 4 g p a ，v 二= o 2 。数值计算利用通用有限元程序a n s y s 图中的模型被划 分成1 4 1 9 7 个四面体单元，共有2 2 1 9 7 个节点，最小的单元尺寸为5 l o “口，计算 时间大约为3 0 分钟。在单元划分时，由于无法预先按一定的比例在立方体内的某 一方向上布置有关节点，因此数值分析所需的位移分量只能通过插值得到。由第 四章中的特殊有限元法可得到，界面端有三重实应力奇异性，奇异性次数分别为 = o 1 6 4 4 ，也= o4 7 0 7 ，如= o 5 7 0 7 。参照三重应力奇异性的应力解法，用位 移分量代替( 2 2 6 ) 式中应力分量，相应的用p - _ 1 代替( 2 2 5 ) 式中的七，可得 鸣。+ 6 t2 + 以，+ j = o ( 2 3 4 ) 式中 o = “f 4 0 p “；：“p 甜p ) 十“p 0 j ! ) “p 一“p “p ) 十“ 0 p “2 ) 一“f 2 ) “竽) 6 = “_ 0 9 砷一“p k ) + “4 ，0 f ”j ”一0 k p ) + 毋0 p k 一“p “ ) c = “f ”0 1 ”“一“5 1 ，“ 2 ，) + “，0 掣“! ”一”f 2 ，“，) + ”p 0 f 2 k 2 一“f 1 ”i 2 ) d = “燃”一“蹦1 1 ) + “卜“燃2 ) + “婚一“m ! 1 ) ) ( 2 - 3 5 ) 由于是三重实应力奇异性，求解方程( 2 3 4 ) ，可以得到三个实根 = p 2 ，一，两边 取对数得 五j = l + l o g 七，1 0 9 p ( 2 3 6 ) 将插值得到的位移分量代入( 2 3 5 ) 式求出t ，然后再代入( 2 3 6 ) 式，就可以 给出如图2 1 7 所示的双对数奇异性分布，利用最小二乘法可求得奇异性次数分别 为 = o 1 5 6 2 ，【：= 0 4 6 7 3 ，丑= 0 5 8 2 1 。与通过特殊有限元法得到的数值结果相 比较，相对误筹分别为5 0 ，0 7 和2o 。 2 1 浙江大学博士学位论文( 2 0 0 2 ) 一 九= o1 5 6 2 k = o4 6 7 3 一 丸一= o5 8 2 1 l o g ( 日) 3 2 图2 1 6 三维界面端奇异性分析模型图21 7 奇异性分布 2 3 5 注意事项 为了精确地求取多重应力奇异性次数，对于平面问题，要求在奇异点附近的 单元划分时按等比级数布置好有关节点；而对于三维问题，可由从界面端出发的 某一方向上的应力和位移分量进行插值，来得到进行数值分析所需的该方向上按 等比级数排列的点上的应力和位移分量。在利用插值方法求取应力分量时，要注 意相邻两点不能在同一个单元内，因为同一单元内的应力分布已由形函数所决定， 难以求取多重奇异性。在利用有限元法计算接合材料问题时，注意到界面上应力 的不准确性，避免沿着界面布点。 2 4 小结 本章提出了一种利用普通的数值分析结果确定多重应力奇异性( 包括振荡奇 异性) 次数的数值分析方法，并通过对双材料界面端模型的计算，验证了其有效性。 得到了以下主要结论： 1 、可以利用通常的数值计算结果来确定多重应力奇异性次数。 2 、用本方法求得的多重应力奇异性次数，与理论解以及由其它的数值方法得到的 结果能够很好地吻合，表明该方法具有较高的精度。 2 2 o 8 6 4 2 0 0 o 0 浙江大学博士学位论文( 2 0 0 2 ) 第三章奇异性分析中的特殊有限元法的应用 3 1引言 由各向异性材料接合而成的异材结构，其独立的材料常数个数的增加极大地 增大了对此类结构界面端奇异应力场研究的难度。因此，现有的理论分析阻“，2 5 m 】 仅仅局限于一些由退化的各向异性材料( 横观各向同性或正交各向异性) 构成的 比较简单的模型。鉴于实际的应力场总是以多重奇异性叠加的形式出现，准确地 确定各阶应力奇异性次数，对于完整地描述奇异点附近的应力场，进而正确地考 虑其评价方法，具有非常大的实用价值。于是，在对各向异性接合材料界面端附 近的奇异应力场进行理论分析的同时，用数值方法求解奇异应力场的应力奇异性 次数的研究也取得了重大进展，除了上节提出的方法外，p a 鼯a u 等人【，9 7 ，5 6 1 提出 了另一种基于y a m a d a 等人啪1 中给出的奇异变换的数值方法，文献【1 5 3 ，m 15 5 1 则提 出了一种基于最小势能原理的求解三维变形问题的特殊有限元法。特殊有限元方 法的基本思路为：从最小势能原理出发，在极坐标系或圆柱坐标系下，通过分部 积分消除r 项，然后利用有限元方法求解奇异性特征方程，得到满足条件 o r e 以) 1 的全部特征值。 本章从最小势能原理出发，在极坐标系下，在通过分部积分消除r 项，从而把 二维问题降阶为一维问题，导出了求解弹性材料平面变形问题和拟平面应变问题 的特殊有限元列式，研究了多相接合材料界面端和界面角点附近奇异性随材料相 数、材料特性以及接合角的变化。数值计算结果表明，本方法具有计算精度高、 速度快等优点。 3 2 弹性材料的平面变形问题 基于y a m a d a 等人m 】提出的奇异变换，p a g e a u 等人1 3 9 】研究了各向异性弹性接 合材料平面变形情况下的应力奇异性。与上述方法不同本节从最小势能原理出 发，通过坐标变换和变量分离消除r 项，给出了一种求解此类问题的特殊有限元法。 3 2 1 基本理论 在极坐标系p ，曰) 下，界面端附近奇异应力场的分析模型如图3 1 ，坐标系的原 2 3 浙江大学博士学位论文( 2 0 0 2 ) 点d 位于扇形的顶点。边界r 2 由从原点出发的两条径向线0 啁和0 口生成，而边界 r i 是半径为的段圆弧，扇形的半径乇为确保圆弧在奇异应力场的支配区域内的 任意常数。为了分析界面端附近的应力奇异性，位移分量可写成如下变量分离的 形式 “= r 。厂p ) ，v = r l g p ) ( 3 一1 ) 式中，“和v 分别是坐标，和臼方向上的位移分量，p ) 和g p ) 是角函数，应力奇 异性次数1 一丑满足o r e ( 1 一 ) 1 。 图3 1 极坐标系下的二维应力奇异性分析模型 相对于楔形体( o ) 整体而言，图3 1 所示的扇形只是它的一部分，也 就是说f 不是真难意义上的自由边界。当楔形体受到一定的位移场作用时，在扇 形体洲b 的边界一上会产生分布力p = ( p ，胁) 。因此，位移场( 3 一1 ) 应该满足扇 彤伟内明半衡乃崔以及边界i j 柙i j 上午h 天明边界条仟。相应的，最小势能原理口j 表示为 面+ r 。函+ 中。却】。一o ，国+ 岛函k k 目= o ( 3 2 ) 式中 刃= 水一抛帜。卜珈枷蚺卜p 。， 这里，中，“和。，的下标表示偏微分( 如“。= 抛归口) ，占表示变分。q 是为了表 2 4 浙江大学博士学位论文( 2 0 0 2 ) 述方便而引入的中问量，其与应变能密度e 之间的关系为中= 旭。由a 纠。毛= 以及几何方程 。= 雾，e 。= ；+ ；嘉，= ；嘉+ 害一； c 。训 5 ”2 石7 。* 2 7 + 7 石万7 r 一2 7 百万+ 万一了 3 4 可得 中。= 盯钾，中_ = r 盯什，o = 盯，日，巾，= 盯柙，中，= r 盯，p ，中= 仃卯 ( 35 ) 将( 3 5 ) 式代八( 3 2 ) 式，得 面= o( 3 6 ) 盯，= p ，口，8 = p 日在边界r 1 上 ( 3 7 ) 由式( 3 7 ) 可知，由于位移场的作用而在边界f 上产生的力的边界条件自动满足。 应力、应变分量可以简记为 盯i = 仃，盯2 = 盯御，o - = 盯，口，】= 什，岛= 占鲫，毛= ，口 ( 38 ) 将位移分量( 3 1 ) 代入几何方程( 3 4 ) 得 = r 。1 仁+ 五) ，( f = 1 ，2 ，3 ) ( 3 9 ) 式中 互= 0 ，乏= 厂+ g 口，弓= 厶一g ，皇= ，岛= o ，善3 = g ( 3i o ) 这里，下标臼表示求偏微分。将应变( 3 9 ) 代入本构方程得 盯，= ，。1 e + 五童) ，( j = l ，2 ，3 ) ( 3 一1 1 ) 式中 石。= c ，| i t ，a 。= cq ij ( 3 1 2 ) 这里，重复的下标代表求和，c 。为与极坐标系对应的材料常数。将( 3 一1 1 ) 依次 代八( 3 5 ) 、( 3 6 ) 得 面= p 。渤= o ( 3 们) 式中 2 5 p酝 旧 抬 嗽嗡嘲璁硝时 f 一 航 p b 蠊肌 + 一q 妒 一 妒黔 矛 五 面 浙江大学博士学位论文( 2 0 0 2 ) 对于所有的驴和国，有硒= o 。 3 2 2 公式的离散化 在进行网格划分时，采用三节点一维等参数二次单元，单元的每个节点有两 个自由度。单元的形函数可用自然坐标表示为 = 一；( ，一7 b ，= 一矿，= 丢( + 叩b ( 3 一，s ) 应用等参数单元的概念，可得 口= 2 目l + 2 目2 + j v 3 岛 式中，伊是节点坐标。单元内部位移的角函数为 ，p ) = f 置，g p ) = g 。x ，( f = l 6 ) ( 3 一1 7 ) 式中，( x l = “l ，爿2 = v ，= “：，。= v 2 ，爿；= 地，片。= v 3 ) ，“，和v ，是节点位移。 函数f 1 ，g 。与形函数m 的关系为 f “= ，g 。= 。，( f = 1 ，2 ，3 ) 将( 3 1 7 ) 依次代入( 3 一l o ) 、( 3 1 2 ) 可得 西j = 弓。i xr ，67 = u x j 式中 瓦= c ，：( f - + g ；) + c ，( 略一g ，) ，雪。= c 。f ，+ c ，g ， 将( 3 一1 7 ) 和( 3 一1 9 ) 代入( 3 一1 4 ) ，然后在单元q f 。1 上积分，可得 面= 硝，+ 列1 ) + 丑2 m ? 弦， 式中 曙= ：j f 七3 i f ；一i ，j g 。十：j g 小d 8 力= i 始：，一气妒。+ 毫，石一恒，十寞，b 十毫，g ；) d 占 砖) | 一鹏，+ 毫，g 口 于早，可以得到以矩阵形式表示的特征方程 ( 3 1 8 ) ( 3 一】9 ) ( 3 2 0 ) ( 3 2 1 ) ( 3 2 2 ) 2 6 浙江大学博士学位论文( 2 0 0 2 ) ( k 十加十名m 皿= o ( 3 2 3 ) 式中，x 是整体自由度矢量，k 、d 和m 是v 的矩阵，是总自由度数。由 于矩阵k 是奇异的【捌，特征方程( 3 2 3 ) 可以唯一地转换成如下形式 a 姜 = 一0 ，k l 一。 萎 = a 姜) c 。一z a ， 由特征方程组( ： - 2 4 ) ，可得到满足条件o r e “ 1 ) 的全部特征值五。 对于各向异陛材料，由于材料常数是角变量的函数，在数值积分时，必须求 出每一个高斯积分点上的材料常数矩阵，详细的推导过程见附录。 3 2 3 程序的开发 首先读入节点和单元信息，形成单元刚度矩阵，然后按照对号入座的办法把 单元刚度矩阵组装成整体刚度矩阵。单元刚度矩阵的形成主要通过以下几个步骤： l 、采用四点高斯积分，计算与位移角变量相关的列式。 2 、计算每个高斯积分点上的材料常数矩阵。 3 、计算每个高斯积分点上的单元刚度矩阵，然后加权相加，给出单元刚度矩 阵m 驴、卯和t p 。 在计算单元刚度矩阵的过程中关键是确定每一个高斯积分点上与角度有关的材 料常数矩阵。整体刚度矩阵形成以后，需要引入位移边界条件。如果己知位移 “，= 0 ，可分别将整体刚度矩阵k 、m 和d 主对角线上的元素女，m 。和一，换成1 ， 相应的把第f 行和第f 列上其余元素换成o 。然后，将三个整体刚度矩阵合成为个 非对称实矩阵a 。由于矩阵a 是非对称实矩阵，它的特征值可能对矩阵元素的微 小变化很敏感，矩阵本身也可能是缺损的，因此通常采用的求特征值的数值方法 已经不能适用，本文采用一种寻找由非对称实矩阵得到的实海森伯格( h e s s e n b e 职) 矩阵的全部特征值的鲫算法。 3 2 4 数值结果和讨论 l 、与界面垂直相交的裂纹模型 2 7 计i l 大学博士学位论乏( 2 2 奇异性数值分析模型如图32 所示裂纹与界面肼垂直，材利常数吣表3l ， 材料主轴的方向由图32 中的坐标系( 1 2 ，3 1 给出，降阶后得到的维线性钶! 域如吲 33 所示，线段的端点爿和，j 分别表示裂纹表面f 朋和f ) ，内点疗和( 分别裹示界 面f 旭和似1 。线段爿b 和( y ) 表示材料i 所占的领域，而线段疗( 表示材芈12 所占 的领域。 3 叠卜? ，) 一 h 料l 图32 与界面母直相交的裂纹模型 q l 疗 q 2 f q ； ) 口= 0 秽= ，曰= y + 疗 目= 2 丌 图33 一维线性领域 存数值计算过程中，采用四点高斯积分。每种材料所占的线性领域各做分成6 1 、大小相等的单，庄一个单元内只能囱种材 ：i 庄平面应力状态 j 子种材 料组合情况下的汁算结果如表32 所示与已有的理论解相比较，在较籼的刚恪划 分下，仍能达剑很高的精度， 表3l各向异性材料的l 催材 ：l 常数 材料编号 i i 。l ；n m ：、e ，( ( 删脚、)( j ，、( ( 删伽、) ，、 l5 51 6l7 06 548 3u 【j 3 6 21 3 44 53 l0 32 41 306 5 0 3l5 47 71s s8 35 q6 803 0 0 2 8 浙江大学博士学位论文( 2 0 0 2 ) 表3 - 2 数值计算结果与已有理论解的比较 材料编号应力奇异性特征值 q i ，q 3 q 、 d e l a l e e r d o g a n c h e n h a r a d a ( 19 9 7 ) 本文的结果 ( 1 9 7 9 )t 五 l2o 4 4 9 5 20 4 4 9 4 8 o 6 2 3 6 80 4 4 4 9 5o 6 2 0 6 7 20 3 4 3 0 l03 4 3 0 10 4 6 7 6 0 0 3 4 3 6 20 4 6 8 8 0 2 105 7 7 4 2 05 7 7 4 4o 3 5 3 2 7 0 5 8 4 0 80 - 3 6 0 8 8 2o 6 3 0 8 906 3 0 8 8o 5 4 4 1 7 o6 3 1 1 3 05 4 4 6 7 为了验证单j 二形状畸变对计算精度的影响，采用不同的单元网格划分，材料l 中单元的大小仅仅是材料2 中单元大小的十分之一对于表3 2 中的第一种组合， 数值计算得到的特征值为 = o 4 4 9 9 4 和t = 0 6 2 4 2 8 ，与表3 2 中已有的结果相比 较，相对误差分别为1 1 2 和0 5 8 。由此可见，单元形状畸变对计算精度影响非 常小，因此没有必要应用大小相等的单元。 2 、三相接合材料模型 三相接合材料界面端模型如图3 4 ，模型中的材料1 和2 是两种泊松比相同而 弹性模量和剪切模量相差一个数量级的正交各向异性材料，其工程材料常数为： 强度较大的材料 1 = 1 o l o p a ，e ! = 2 e i ，= 3 e i ，g 1 2 = 4 巨，r i2 = p 13 = 屹2 = o 3 强度较小的材料 1 = 1 o 1 0 8 p a ，e 2 = 2 e 1 ，毛= 3 e l ，g 1 2 = 4 e l ，v 1 2 = v 1 1 = v 3 2 = 0 - 3 材料的弹性主轴如坐标( 1 ，2 ，3 ) 所示。角度口的变化范围为：o o 口蔓9 0 0 ，模型降阶 后得到的一维线性领域被划分为1 2 个单元，采用四点高斯积分。当材料l 的强度 较大时，特征值刖五，u = 1 ，2 ) 表示：反之，特征值用a ，o = 3 ，4 ) 表示。 在平面应变状态f 界面端0 附近奇异性次数随角度口的变化曲线如图3 5 所 示，由奇异性变化曲线可以发现：当材料1 的强度大于材料2 时，界面端附近的 应力奇异性次数1 一五，( ，= l ，2 ) 随着角度口的增大而递增，在口= 2 6 0 附近，奇异性 2 9 浙江大学博十学位论文( 2 0 0 2 ) 由一个增加为两个，当口= 4 5 0 时较强的奇异性出现一个极点，当口：9 0 0 时奇异性 最强。当材料1 的强度小于材料2 时，在角度口 1 0 0 的情况下界面端不产生奇异 性，在a = 1 0 0 附近出现第一重奇异性，在口= 5 0 0 附近出现第二重奇异性，当口：9 0 0 时奇异性最强。对于任何角度口，白异性次数卜丑，o = 1 ，2 ) 总比1 五，( j = 3 ，4 ) 要强 得多，这可以理解为相对于强度较小的材料，另一种材料可看作是它的刚性约束。 jiy 材料2 0 x 槲筋受荆 图3 4 三相接合材料界面端模型 图35 奇异性次数l - a 随接合角口的变化曲线 3 0 浙江丈学博士学位论文( 2 0 0 2 ) 3 、四相接合材料模型 四相接合材料界面端和界面角点模型如图3 6 ，材料常数由上节给出，弹性主 轴如图中的坐标( 1 ，2 ，3 ) 所示。角度口的变化范围为：o o s 口1 8 0 0 ，降阶后得到的 一维线性领域被划分为1 2 个单元，采用4 点高斯积分。 千2jiy 1 一 ) 、 稿料2 l 舀厂 材料1 r 7 aj 图3 6 四相接合材料界面角点和界面端模型 界面端d 附近奇异性次数1 一五随角度口的变化曲线如图37 所示，当材料1 的强度小于材料2 的强度时，奇异性次数用l 一 ，o = 1 ，2 ，3 ) 表示，反之，奇异性次 数用1 一a ，u = 4 ，5 ，6 ) 表示。由奇异性变化曲线可以发现：两种情况f 的奇异性关于 d = 9 0 0 对称，最多有三个奇异性，最强的应力奇异性l 一工= o 7 4 4 2 分别出现在 瑾= 1 2 5 0 和口= 5 5 0 附近。当口= o o 或d = 1 8 0 0 时，图中的模型退化为同种材料中的 裂纹模型，裂纹尖端的奇异性次数与材料常数无关，两重奇异性次数皆为o 5 。 界面角点d 附近奇异性次数l 一丑随角度口的变化曲线如图3 8 所示，当材料1 的强度小于材料2 的强度时，奇异性次数用1 五，( j = l ，2 ) 表示，反之，奇异性次数 用卜兄，( u ，= 3 ，4 ) 表示。由奇异性变化曲线可以发现：两种情况下的奇异性关于 口= 9 0 0 对称，最多有两个奇异性，最强的应力奇异性l 一五= o 6 1 1 1 分别出现在 3 l 浙江大学博士学位论文( 2 0 0 2 ) a = 1 2 5 0 和口= 5 5 0 附近，这一点与界面端附近的奇异性恰好相同。值得指出的是， 当口相同时，界面端附近最强的奇异性都比界面角点附近最强的奇异性大，奇异 性的个数也较多。 1 一九 图3 7 界面端附近的奇异性次数1 一 随接合角口的变化曲线 l 一九 图3 8 界面角点附近的奇异性次数l a 随接合角口的变化曲线 3 2 浙江大学博士学位论文( 2 0 0 2 ) 3 3 弹性材料的拟平面应变问题 对于图39 所示的模型作用在其上的载荷仅仅是坐标x 和y 的函数，而与坐 标z 无关。相对了二x 和y 方向，z 向的尺寸足够大，从而确保应力和位移也与坐标z 无关。由于材料是各向异性的，( x ，y l 平面不再是一个材料对称面。尽管是一个二 维问题，但在变形后二维平面将成为一个三维曲面，应力状态是三维的。这类问 题可以称做拟平面应变( g e n e r a l i z e dp l a j l es t r a i n ) 问题。基于y a m a d a m ”提出的奇 异变换，p a g e a u 和b i g g e r s 一”研究了各向异性弹性接台材料界面端和界面角点拟平 面应变情况下的应力奇异性。与上述方法不同，我们从最小势能原理出发，通过 坐标变换和变量分离消除，项，给出了一种求解此类问题的高效的特殊有限元法。 图3 9 拟平面应变问题模型图3 1 0 圆柱坐标系下的奇异性分析模型 3 3 1 基本理论 奇异性分析模型如图3 1 0 所示，在圆柱坐标系( r ，日，z ) 下，位移分量可写成如 下变量分离的形式 “= ，j p ) ，v = ，。g p ) ，w = r 2 p ) ( 3 2 5 ) 式中，“、v 和w 分别是坐标r 、口和z 方向上的位移分量，p ) 、g p ) 和 p ) 是 角函数，应力奇异性次数i 五满足o r c ( 1 一五) l 。 3 3 - 浙江大学博士学位论文( 2 0 0 2 ) 类似与弹性材料平面变形问题的分析，可以认为在边界f 上作用有分布力 p = b ? ，p ：p ? ) ，位移场( 3 2 5 ) 应该满足扇形体内的平衡方程以及边界r l 和上 相夭嗣】卫界条忏。干日比酬，最，j 、劈日e 原埋口j 衣不刀 办+ r 。面+ 中。西+ 巾。丽】，扫? 函+ 硝函十p ? 辄k k 口= o 式中 面= 小一鼢帜坼+ 卜珈枷叩。卜 + 心一抄帖吣帆卜 式中，巾= 饱( e 是应变能密度) 。由钮a ，= 盯。以及几何方程 ( 3 2 6 ) ( 3 2 7 ) 知= 雾，= ；+ ；嘉，= 吉嚣+ 害一；，靠= 警，= ；荔c 。瑚，知5 石2 了+ 7 丽，一2 7 历+ 石一7 ，一2 i 2 ；历叫w 可得 0 。= 仃鲫，m = r 仃，由如= 盯，口，中p = 一盯坩，中= r 盯，口 m = 盯舳，o 。= o ，o 叶= r 盯。，中= 盯 将( 3 2 9 ) 式带入( 3 2 6 ) 式，得 厕= o 仃，= p ? ，盯，。= p ：，盯。= p ? 在边界l 上 ( 3 2 9 ) ( 3 3 0 ) ( 33 1 ) 由式( 3 3 1 ) 可知，由于位移场的作用在边界f 上产生的力的边界条件自动满足。 应力应变分量可以简记为 盯l = 盯什，仃2 = 盯鲫盯3 = 盯出盯4 = 盯口，盯s = 盯r 日 q = 占。，2 = ，f 3 = ，t 4 = y 。，毛= ，日 ( 3 3 2 ) 将位移分量( 3 2 5 ) 代入几何方程( 3 2 8 ) 得 式中 = ，“1 仁+ 丑皇) ，( f = 1 5 ) i l = q ti ! = j + g9 ，i l = h 9 ，i = o ，i 5 = j 8 一g ( 3 3 3 ) 3 4 浙江大学博七学位论文【2 0 0 2 ) 善l = - 厂，0 2 = o ，言3 = o ，善4 = ，言5 = g 将应变( 3 3 3 ) 代入本构方程得 式中 q = r 。1 ( e + 五号) ，( f = 1 5 ) 否? = c i ij ，6 = c 。； ( 3 3 4 ) ( 3 3 5 ) ( 3 3 6 ) 这里，重复的一f 标表不求和，c 。为与圆柱坐标系对应的材料常数。将( 3 3 5 ) 依 次代入( 3 2 9 ) 、( 3 3 0 ) 得 办= i r 2 “渤= o( 3 哪) 式中 撕= 一f p 形+ 也国+ 号。踊) d 目 + f 陋：一瓦阿+ 子，诉一眩+ 子，壤+ 子z 哦一瓦劢+ 方，】d 口 ( 3 3 8 ) + f 慨酽+ 瓦矾一元老+ 五露。+ 五劢。) d 口 对于所有的可，国和丽，有撕= o 。 3 3 2 公式的离散化 在进行网格划分时，采用三节点维等参数二次单元，单元的每个节点有三 个自由度。单元内部位移的角函数为 p ) = ，。置，g p ) = g z ，厅p ) = z ，( = 1 9 ) ( 3 3 9 ) 式中，( x 。= “，2 = v 。，扎= w 1 ，。= “2 ，= “，一) ，“。，v ，和w 是节点位移。 函数f ，g 。，h 与形函数旷的关系为 ，“! = ，g ”1 = ，h “= _ v ，( i = l ，2 ，3 ) ( 3 4 0 ) 将( 3 3 9 ) 代入( 3 3 4 ) ，可得 i ? = e ，x 。，? = e 。xj ( 3 4 、 式中 置1 i = o ，重! 。= fj + g ；t 毛3 i = h ；t 毫4j = o ，盖s i = f ；一g i 3 5 浙江大学博士学位论文( 2 0 0 2 ) e ，= f te 2 ，20 ，毛，= 0 ，e ，= h7 ，巨，= g ( 3 4 2 ) 将( 34 1 ) 代入t3 3 6 ) ，得 互= 写，声= 毫， ( 3 4 3 ) 式中 sq = c t e bts u = c e b ( 3 4 4 ) 将( 33 9 ) 和( 一 _ 4 3 ) 代入( 3 3 8 ) ，然后在单元q 上积分t 可得 痧= 料，+ 列”) + 五2 m p 汪， ( 34 j ) 式中 m ? = 一l 幢，f 。+ 寞，g + 文，7 ) d 口 d p = f 呤：，墨，p 。+ ，g 一慨，+ ，) g + j ：，q 一只，h 。+ 毫，蟛) d 毋 t = f 恒，f + 瓦，巧一豆，g + 夏，g ；十墨，h ；) d 口 ( 34 6 ) 于是，可以得到( 3 2 4 ) 所示形式的特征方程组。圆柱坐标系下高斯积分点上的 材料矩阵详细的推导过程见附录。 3 3 3 数值结粜和讨论 1 、与界面垂直相交的裂纹模型 纤维加强复合