（机械设计及理论专业论文）斜齿轮接触问题的形状优化研究.pdf

太原理工大学硕士学位论文 摘要 本文以边界元法作为应力分析工具，对三维弹性接触 问题的形状优化问题进行了理论上的探讨，并对斜齿轮接 触问题进行了实际分析。作者编制了斜齿轮前处理程序， 应用边界元应力分析程序3 d c b e m 分析了沿齿宽方向的 轮齿接触应力的分布。选取约束变尺度法作为优化算法， 对斜齿轮进行了齿廓优化，得到了满意的结果。 关健词：边界元形状优化斜齿轮接触问题 第2 页 太原理t ：大学硕士学位论文 a b s t r a c t t h i sp a p e rd i s c u s s e st h eo p t i m i z a t i o np r o b l e mo f 3 一de l a s t i cc o n t a c to b j e c t si nt h e o r yu s i n gb e m ，a n d a n a l y z et h er e a lc o n t a c tp r o b l e mo fh e l i c a lg e a r s t h e a u t h o rd e v e l o p e dt h ep r o g r a mw h i c hi sn e c e s s a r yf o r t h es t r e s s e s a n a l y s e s a n a l y z e t h es t r e s s e s d i s t r i b u t i n go fh e l i c a lg e a ri nl o n gd i r e c t i o nu s i n g 3 d c b e m o p t i m i z et h es h a p eo fd e n t i f o r mu s i n g c v ma st h eo p t i m i z ea l g o r i t h m t h er e s u l tw e r e s a t i s f i e d k e y w o r d s ：b o u n d a r ye l e m e n tm e t h o d s h a p eo p t im iz a t i o n h e l i c a lg e a r c o n t a c tp r o b le m 太顾理工大学硕士学位论文 第一章绪论 弹性力学中研究两接触物体受力后产生局部应力和 变形的问题统称为接触问题。工程结构中，系统往往可分 成几个永久或非永久联接在一起的部分，这些部分之问的 力是靠它们之间的接触挤压，甚至足冲击柬传递的，因此 接触问题在工程实际中经常碰到，例如：齿轮的齿问相互 啮合接触、凸轮与挺杆的接触、传动轴与轴座的配合接触、 轴承的滚珠或滚针与滚道的接触等。这些零件的寿命大多 与其接触面上的挤压和接触应力的大小和分布有关因 此，确定接触区域的大小及接触压力的分布规律就成为工 程技术领域一个重要的课题。国内外在这一方面有很多文 献，我室多年来一直从事这方面的研究，取得了很大的成 绩。 工程中的接触问题是各种各样的，主要有：弹性物体 的接触、塑性物体的接触、粘弹性物体的接触、可变形固 体同液体之间的接触。本文所研究的接触问题足弹性物体 的接触。接触问题有两种求解方法：解析法神1 数值法。接 触问题的解析法又称经典接触力学。经典接触力学分析接 触问题只局限于一些几何形状比较规则的物体，应用范围 非常有限工程中遇到的实际接触问题，其结构形状、边 界形状及接触状态都十分复杂，一般无法求得解析解。接 触问题的数值解又称为非经典接触力学。近几十年来，随 着应用数学、应用力学理论和数值计算以及计算机技术的 迅猛发展，特别是有限元法、边界元法等数值方法的出现， 使得非经典力学得到越来越广泛的应用，它为我们分析工 程实际中复杂的接触问题提供了强有力的工具，接触问题 第4 页 太原理工大学硕士学位论文 的研究取得了很大的发展。数值解法主要有：有限元法、 边界元法、有限差分法以及与数值方法相配合的各种变分 法、实变函数法、泛函分析法等。 目前，人们分析接触问题的工具主要有有限元法 ( f e m ) 和边界元法( b e m ) 。有限元法产生于5 0 年代初， 它把差分法的离散改造成更为灵活的有限单元离散，把早 兹法的试函数近似换成插值函数近似，以弹性力学变分原 理作为推导的根据，尤其是它充分利用了电子汁算机的计 算能力从而开拓了现代数值方法的广阔领域。特别是有限 元简单直观，易于掌握而且使用范围广，计算效率高， 数值精度也较好。人们已用它来求解各种力学和非力学问 题、线性和非线性问题，均取得好的成效。目前已经有许 多相当完善的程序和程序包，诸如s y n s a s 、a d i n a 、 s a p 9 l 和s a p 9 5 等等，已广泛应用到结构件的应力分析和 结构优化设计中。尽管有限元法所取得的成就与日俱增， 但有限元法还不是十全十美的。改进有限元法的努力一直 在进行着，但有限元法的某些不足是无法克服的。例如有 限元法需全域离散，导致问题的自由度和原始信息量大： 对无限域只能人为的取成有限域：有限元法的离散技术本 身也存在着缺陷，它把本来是连续的介质用仅在节点处连 接的有限单元的集合来模拟，这样不仅带进了离散的误 差，而且在单元之间连续的要求较高时，有限单元的构造 也很困难；对有限元法的精度和可靠性也常常会提出疑 问，因为对同一问题采用不同的程序计算时可能会得出不 同的结果。 边界元法，实际上是指边界积分法及其数值解法。边 界元法可分为间接法和直接法，现在工程界常用的是直接 法，即用加权余量法推导边界积分方程。得到的积分表达 第5 负 太原理工大学硕士学位论文 式是区域内物理量与边界上物理量之间的关系式，积分方 程中的未知函数是所求物理量在边界上的值。加权余量法 的基本思想是把近似解和真解在区域和边界上的误差与 某个选定的权函数构成数量积，使得误差在区域和边界上 按一定方式分配从而使误差趋于极小。边界元法是选取 在区域内满足控制方程，且近似地满足边界条件的解作为 权函数，或者选取任何不带任何边界条件的奇异解作为权 函数，也就是以基本解作为权函数。边界元法同时又包含 有限元法的思想，它把有限元法按求解域划分单元离散的 概念移植到边界积分方程方法中。但边界元法不是有限元 法的改进或发展，边界元法与有限元法存在着质的差异。 有限元法要在整个求解域上进行离散，边界元法只在求解 域的边界上离散；有限元法是完全的或全域的数值方法 而边界元法在域内采用了弹性力学基本解和s o m i g l i a n a 积 分，数值计算只在边界上进行，它属于半解析半数值方法。 不难看出边界元法具有有限元法所没有的优点。 有限元法的不足看来只能用边界元法来弥补。边界元 法仅在边界上离散，使数值计算的维数降低一维，从而减 少了问题的自由度和原始信息量。边界元法采用无限域的 基本解，用边界元法求解无限域问题可称是天衣无缝。边 界元法的离散误差只产生在边界上，边界元法中部分采用 数值法，部分采用解析式，它的准确度和可靠性已公认是 高于有限元法。对于接触问题，确定接触区域的大小和接 触压力的分布规律的问题，本身就只局限于问题的边界 上，故采用边界单元法比有限元法更为有效。 形状优化就是合理地设计机器构件的形状以提高构件 的性能和功能，它是设计过程中较高层次的问题。形状优 化设计虽然是近年来才兴起的研究领域但已有不少成 第6 顶 太原理j t ：大学硕士学位论文 果。在结构分析和机械设计的应用中有合理选择承载件 的形状，以使构件某方面性能最优如重量最轻、最大应 力极小等等。在接触问题的应用中，有接触应力极小化的 轮廓设计、支承轮廓设计等等。边界元法和有限元法是建 立形状优化数学模型的有力工具，与有限元法相比，边界 元法在形状优化中具有突出的优点，如只需离散问题的边 界，降低问题的维数，计算精度高网格划分简单，单元 数少，占用内存少，计算速度快等。因而。边界元法在形 状优化结构分析中具有更大的优势。但是边界元法目前不 象有限元法那样拥有完善的前后处理功能，因此用边界 元法进行形状优化仍是一个处于不断发展中的研究领域。 近十多年来，我教研室在该领域进行了大量的研究活动， 发表了一批论文。文献【3 l 】对受扭回转体用边界元法进行 形状优化作了研究：文献【3 2 以边界元法为应力分析工具， 以三次样条函数描述待优化边界对阶梯轴轴肩进行了形 状优化，并用解析法推导了应力敏度分析公式：文献【3 3 1 建立了用边界元法进行形状优化的前、后处理及优化过程 智能化处理的模式；文献( 3 4 】伴随变量法与直接法相结合 的敏度分析的混和方法，该方法提高了敏度分析的效率 解决了三维形状优化机时过多的一大困难：文献【3 5 】提出 了一种新的敏度分析方法一变分法，减少了优化迭代次 数，提高了计算速度。文献【2 8 】对二维弹性接触问题的形 状优化，进行分析和整个系统设想以及算例分析。文献【2 7 】 对三维接触问题形体的形状优化作了理论上的探讨，并自 行开发了使用于简单形体的边界元应力分析程序。文献3 6 1 应用边界元法对直齿圆柱齿轮的应力进行计算，并且在应 力分析的基础上进行了形状优化。这些文献对用边界元法 进行形状优化的理论研究和工程应用作出了一定的贡献。 第7 页 太原理工大学硕士学位论文 尽管人们在研究用边界元法进行形状优化方面取得 了不少成果，但是在应用边界元法对斜齿轮进行形状优化 方面的文献很少本文在参阅大量有限元斜齿轮网格划分 的基础上编制了斜齿轮前处理程序，运用文献r 2 7 的应 力分析程序对斜齿轮接触应力进行了计算，并对斜齿轮齿 廓进行了优化，得到了满意的结果。 作者主要做了咀下工作： 1 、详细阐述了三维弹性力学的边界元法和形状优化 的算法原理。 2 、 详细阐述了斜齿轮建模所需的基础知识。 3 、 编制了斜齿轮前处理程序，只需输入齿轮的基本参 数、材料系数等参数便可求得各个节点坐标和编号，自动 进行网格划分。 4 、 修改了应力分析程序，大大减少了原始数据的输入 量，提高了边界元分析应力的效率。 5 、对斜齿轮进行了齿廓形状优化，得到了满意的结 果。 第8 页 太原理工大学硕士学位论文 第二章带摩擦三维弹性闩曩的边界元法 2 1 边界积分基本方程 1 8 4 8 年，开尔文求解了集中力作用于无限大物体内的 点上的问题，得出了弹性力学基本解k e l v i n 解。贝蒂于 l8 7 2 年提出了互等定理。s o m i g l i a n a 把b e t t i 互等定理和 k e l v i n 解相结合，对待求解和k e l v i n 解采用了b e t t i 互等 定理得到了用边界位移和面力表示内点位移的边界积分公 式，即：s o m i g l i a n a 恒等式。弹性力学边界积分方程的推 导和s o m i g l i a n a 恒等式的推导相似，两者的不同只在于现 在需将源点移到边界上，得到的弹性力学边界积分方程为： “：+ f p 以d f = i ”a p 。d f + ”：以拉 ( 2 1 ) 其中，= l ，2 ，3 ：蜥与p 。为边界r 上的位移与面力张量；“m 和肌为位移和面力的基本解- - k e l v i n 解。k e l v i n 解的基 本形式为： “池g ) = 面石1 面阶4 v ) 毛+ 燮 忡，= 高御功，c 盟掣 【( i 一2 y ) d 。+ ! ! ! ：二! ! ：! ；! ! ：_ = 翌1 1 ! ! ：望【( i 一2 l ，) d 。+ 二l i 土一1 上二二 口v v 其中，为场点p 到源点q 的距离，v 为泊松比，“：，p ：的 第9 弧 太腺理工大学碗士学位论文 第一个下标i 表示p 点在x i 坐标方向的位移和面力分量，第 二个下标j 表示作用于源点q 的单位集中力指向x j 坐标方 向 由于f c 。万d r = c m ：( 占。为d i r a c 函数) ，因此式( 2 i ) 改写 成： f p 订= i “：仇d r + b ； ( 2 2 ) 其中五= 五十c 。万；4 为( 2 1 ) 式中的体积力积分项。 2 2 三维弹性接麓向息的边界积分方程 2 2 1 接向曩的基本假定 1 接触体之间无刚体位移 这一假设是为了消除由于刚体位移而导致的总体矩 阵的奇异性，以简化问题的处理，对于必要的边界条件总 是可以办到的。 2 接触面光滑。摩擦力服从库仑定律 光滑接触面的假设，给问题在数学f = 的处理带来了方 便。对于金属材料这一假设是合理的：用库仑定律来确定 摩擦力较为符合实际情况。 3 接触体为线性各向同性材料，且为小变形 这条假设保证了线弹性及小变形理论的使用。在不太 大的外载作用下。这条假设总是成立的。 4 略去接触面问杂质( 如油、粉尘等) 的影响 接触面间的问距一般说来是很小的。即使有少量的 油、粉尘侵入，对接触面上力学性质的影响也是很小的。 在该假设的前提下，可以略去由杂质而导致的耗散性质。 第1 0 页 太原理工大学硕士学位论文 2 2 2 增量迭代法 对于接触问题的计算，根据接触过程的特征，要做如 下三方面的考虑。 第一，几何非线性问题。这主要表现在接触面积的大 小是随外载荷变化的。在数值分析中要解决的实际问题是 弹性体有关边界从接触到分离和从分离到接触的变化。 第二，接触状态问题。这时主要是指滑动接触和粘着 接触，在数值分析中要解决的是滑动与粘着的互相变化和 滑动方向这两个实际问题。 第三，不可逆问题。 这三个问题中前两个可用迭代计算进行处理，而第三 个问题实际上要体现出加载过程，因此可用增量法来解 决。综合起来就是增量迭代法。 不可逆的主要特征是在某外载荷状态以后，摩擦力 方向有向相反方向变化的趋势。而这一特征通常是在外载 荷逐渐增加和逐渐减小的转变处产生。为了既能反映不可 逆性这一特征，又能节省计算时间，当计算接触体在某 载荷下的解时，可按如下方法选择增量的大小。 l 、若某载荷从初始值p o 以逐渐增大的方式变到 p l 时，可取一次增量，其大小为舯= 只一异；当然也可以 取多次增量，但算倒表明两种情况计算结果完全一致， 与增量大小的选择并无直接关系。 2 、若某一载荷p 从初始值p o 以逐渐增大( 或减小) 的方式变到p i ，随后又以逐渐减小( 或增大) 的方式变 到p 2 ，这时可取两次增量，其大小分别为尸= e 一只和 a p 2 = 只一鼻。 第l i 页 太原理l 人学硕士学位论文 当有多个载荷作用于相互作用的弹性体时，增量的选 择不仅要按上述原则，还要看它们之阃的相互关系。对于 载荷变化较复杂的情况，按上述原则类推，可能要取两个 以上的增量。 按上述方法选择增量时，由于增量个数少。因此，在 其他条件相同时，可以花费较少的计算时问。当接触体间 无摩擦时，可取一次性增量。且仅存在几何非线性，因此 只需对接触区的面积大小作迭代计算。当考虑接触体问摩 擦，但无加载说明时，增量迭代和全量迭代得到相同的结 果。 总的量用和来计算即位移和面力为： u j “= “? + “f = a u ? t = l ( 2 3 ) = 卸? + p = e 卸? ( 2 4 ) k = l 2 2 3接触条件 两弹性体接触时，其边界可分为非接触区和可能接触 区，可能接触区比实际接触区大，通常由粘着区和滑动区 组成。不同的接触区域有不同的边界接触条件。 我们以两个物体相互接触为例来说明各种接触状态 的接触条件，如图2 k a ) 给出两个弹性体q 。和q 8 ，用r 。 和r 8 分别表示两物体的边界。为叙述方便，以下用盯代表 a 和b ，p 为可能的接触边界( 即含部分自由边界) ，的 点q 。和硭的点q 8 构成可能接触的点对。由于接触问题发 生在局部边界上，其接触条件需要以接触面的局部坐标来 表示。因此，需要在接触面f 上建立局部坐标系( 备，喜：，岛) 。 在图2 1 ( b ) 中。坐标系选取如卜： 第1 2 页 太原理工大学硬士学位论文 p 一_ 0 ) r f c 甜 t 0 。 k 。一1 一一譬0 、“ 叩f ，口 芦一i 卜醛 1 日 ( b ) 圈2 1a 和b 物体接状杏圈示 ；：= ：= t ： 芎? = 一 ； 如果采用载荷增量法在小变形的条件下 状态下的接触条件可以表示成： ( 2 5 ) 各种接触 ( 1 ) 分离状态一一相对应的表面处于未接触状态t 并 把该接触面边界记为r ，则接触定解条件为： 幻? = g ? = o ( j = l ，2 ，3 ) 而且满足不等式条件( 判别条件) ，即 一一v 5s v 0 4 3， ( 2 6 a ) ( 2 6 b ) 其中q 、v 表示在局部坐标系下的面力与位移，v 为第打步 载荷增量迭代时接触点对之问在蠡方向的距离，并且 = 陋i 一心川一v ? 一v ；为接触点对之间的初始距离参方 向的分量。 第1 3 页 甜7 当 一 1 太原理工大学硬士学位论文 ( 2 ) 牯着状态一一保持接触。对应表面之间没有相对 位移把该接触面的边界记为聪，则接麓定解条件为： + 缸：= v y 却? 一却? = o + v 8 = i k , o - m 却。一幻? = 0 (27a)2 22 1 2 、 一a v b ，一_ _ _ v ，0 幻；+ 叮；= o331 ，。3 而且满足不等式条件 ( 2 7 b ) ( 3 ) 滑动状态一一保持接触，在接触面之间有相对位 移。把该接触面记为砰，则接触定解条件为： 0 一吣= v ? 。幻? 一幻? = o】，。1 l 幻? = 一g ? 厶g ；一幻：= o 厶叮：= 声：蛔? 幻? + 幻；= o 而且满足不等式条件 ( 2 8a ) 第1 4 页 j 2 ，。：，_-i-f_-_l 列口3知 如卸阻b e k k 瞄瞄 p 肛 p j j = 一2 脚 酗 ，i-ii-f，_-_ 太原理工大学硬士学位论文 垆蛐( 2 8 b ) i 叮；0 阳和胁分别为晶和善z 方向上的摩擦系数，符号“+ ”和“- ” 可根据接触点的滑移方向相反的原则选取。 在一定的外载作用下，每对接触点都必须处于三种接 触状态之一，且满足接触条件中的等式条件。求解时须先 假定接触状态，即根据上述接触等式条件建立方程求解， 结果应满足不等式条件，否则修改接触状态并相应修改接 触等式条件再求解直到在可能接触区内的全部接触点都 满足接触等式条件和不等式条件。因此，可以看出，接触 问题的求解过程是迭代求解过程。 2 2 4 加载的不可逆性 两个接触体之间的压力分布可以是各式各样的而且 一般说来其压力分布是非线性的。但在许多情况下，加载 仍是可逆的此时压力分布与载荷之问呈单值函数关系 载荷一定接触压力也唯一地被确定，即仍具有加载可逆 性。 与此相对照，如果即使载荷相同但加载顺序与方法f 加 载履历) 不同，则其压力分布、应力、变形等均不同，这种 情况称为加载的不可逆过程。 如图2 2 所示结构，我们看到b 物体上s 点的接触压 力，法向接触压力为p 是一常数，水平( 切向) 接触压力q 不仅和q 有关，而且和是否存在p 以及和q 的加载次序 有关。 第1 5 页 阻p 习 r 里一旦些 习。lq 熏h i e q - - 蕊 圈2 2 设想如果只加p 力不加o 力，显然q = o 。 如果只加q ，不加p ，则a 、b 表面之间无约束作用， 则p = o ，q = 0 。 如果先加p 再加q ，然后卸除q ，则其最终载荷状态 和第一种情况一样，但其接触力的状况却不一样，这可从 圈2 2 ( b ) 看出对应接触点的水平接触力q 的变化。 图2 2 ( b ) 表示在p 不变的情况下，q 随q 加载变化而变 化的过程：a b c 段表示q 力从0 加载到q o 过程中q 的变 化：c d e 段表示q 力从q 卸载到0 过程中q 的变化。 从图中可以看出，当q 从0 加至q o ，然后从q 。减至 0 时，接触点处的切向力并没有回到初始值p 2 0 , 而是另一 个值up lo 另外在某一个q 值下，切向有两个可能值， 其中一个对应于加载过程的值p 2 + ，另一个对应于卸载过程 的值p 2 。这种现象就是不可逆性。产生不可逆性的原因 是由于接触面存在摩擦，在往复加载的过程中摩擦力的方 向也要做相应的变化，这就使得接触体的受力与加载有 关。凡是能够出现接触内力不可逆的加载过程可称之为复 杂加载；否则称为简单加载。 既然接触状态的可逆性与摩擦力有关，那么无摩擦接 触自然属于可逆过程：对于有摩擦接触是否可逆，还要看 第1 6 页 太原理工大学硕士学位论文 其具体的载荷情况。对于简单加载的可逆过程可以用增量 解即一次加载求解t 面对于复杂加载的不可逆过程则必须 用增量法求解。才能真实反映接触面切向力的变化状况， 即在接触面迭代计算中对不同性质的载荷按加载的历程 多次麓加于物体上求解。 值得一提的是，在许多接触物体的计算中，采用了把 总的载荷向量分成几级加载的求解方式，一般也把这种方 式称为增量加载，但这和我们讨论的不可逆接触问题的增 量加载是不同的，实际上是不能处理接触问题的不可逆计 算的，但采用此方式的好处在于使载荷分级地施加于物体 上，这样接触点对的状态将平稳变化使求解稳定收敛， 为加以区别把此种加载方式称之为分级加载。在一般的接 触问题计算中可同时采用分级加载和增量加载两种方式。 既可解决接触问题，也可使迭代计算稳定进行。 2 2 5 坐标系和摩攘系数的选取 叩， a ； 聆丽 啦 l 一 球 、f 、立 一一哺 i a 聃 i j | i l 一0 j 。卢 盈2 3在单元第j 节点处的外法线磊圈示 第1 7 页 太原理工大学硕士学位论文 局部坐标系的选取 ( 2 9 ) 如圈2 3 所示，在节点j 处的法向珈按( 2 9 ) 式选取。 但通常坐标系( 协，o r l o ，7 ，钟却：) 不是正交的。当考虑摩擦 力后，必须选取正交坐标系才能使问题变的简单一些。 因此，在图2 1 ( b ) 中，可选取正交坐标系如下： l _ _ a _ 日卜_ 舟f一口 l 季，= ( 叩，一叩，1 j 7 ，一叩，l = 毒， 忙= c 苦+ 篆旃斗篙+ 蕞征i l ； c z 加， 【手? = ( 手：f - - ，a ? ，降：手：i = 一手? 摩擦系数 设在右，蠡平面上总的摩擦系数为u ，如图2 4 所示， 可得到下列方程： g q ：l = ： 御t q 3 ，c 丽o n s 口f i ：= p t t ：c 2 t - ， 如果方程式( 2 8 a ) 和( 2 8 b ) 中的“”号只取“+ ”号， 则有朋= u c o s a ，肫= p s i n a ，并得到下列关系： s i n a = q ：厅可c o s 口= q 抒可 ( 2 1 2 a ) 或者 s i n 口= - - v 2 v ，吒c 。s d = 叫，拇+ y ： ( 2 1 2 b ) 第1 3 页 爵一魄 品一帆 ，m， 茹瓦 x 磊一啊j，j、 = 一以 太朦理工大学坝士字位论文 因此有 _ = 脚：历ip ：= 脚霄可 ( 2 1 3 a ) 或 一= 一：两心= 一一孵可( 2 1 3 b ) 岛 ，i口 卣 v l l 一卜 】 。 v f + v ； 圈2 4 在鲁j 善：平面上的面力与位移关系图示 2 3 6 边界积分方程 如果采用载荷增量法，则( 2 2 ) 式可写成： f p ：缸订= f “：印。订+ a 可 ( 2 1 4 ) 若在接触状态下，需转化为局部坐标，绝对坐标系下 的位移与面力增量与局部坐标系下的位移与面力增量的 关系为： f 缸= 0a v ( 2 1 5 ) 【卸，。吒q 。 其中v t 和q t 是局部坐标系下的位移与面力；靠为嘶与 v 。之间的夹角的余弦值。 石卜崭l 一 a 原理1 人学坝十学何论文 j i i i i i i i i i i i i i i i i ii i i i i i i i i 将力程式( 2 1 4 ) 分别用于a 和b 两物体，并利用 ( 2 15 ) 式，得： 对物体a ： f p ：幽。订+ f 、l q ：舢。d r 2f 。“：印。d r + t v ：砷。打+ 口：( 2 1 6 a ) 对物体b ： f 。p ：l l 。汀+ f l9 ：舢。们2t “：却。d f + fv ：幻。d i g + 曰f ( 21 6 b ) 将接触条件( 2 ，6 ) 和( 2 7 ) 代入( 2 1 6 ) 式，得： 埘物1 奉a ： f i p ：舢汀十fq ：蛳。d 1 1 十lq ：vd f = f ”：卸。d f + iv ：q 。讲1 + f ( 芦，v i 芦! 一! i , ) 6 q ，d f + a b ； ( 2 1 7 a ) 埘物体b ： f 。p ：f f 。汀十卜g j v ? 一口? + q 3 ( a v ? 一v 翔订+ fl q j , a 一 + 9 + g p ? 脚2 ， d f + g ? “i ，a q 2 一v j 知一打十j ( l v i - + 2 v ：w t ) a q ，d 1 1 + 彰 ( 2 1 7 b ) ( 21 7 ) 式即是求解一：维弹性接触问题的边界积分方程的 增量形式。将和r 用单元离散，联立迭代求解，即可求 得边界上的面力、位移及诈确的接触边界r f 。然后，通过 内部点的边界积分公式和边界应力关系式，可分别求得a 、 r 两物体内部点的位移、应力等。 籀2 0 负 太原理工大学硕士学位论文 2 3 敦值技术 2 3 1 积分方程的膏傲 由于直接求解( 2 1 7 ) 式比较困难，所以一般都是利 用数值计算方法求解其近似解。数值逼近法可分为三步： 1 在弹性体的边界上划分网格，网格单元称为边界单 元，单元上的位移和面力取该单元上节点的分段插 值函数。 2 方程( 2 1 7 ) 的离散化的形式应用于边界上的每 个节点，在每一个节点上计算积分。 3 求解上述方程组时，加上己知边界条件( 已知的n 个节点位移和节点力) ，则方程组是一个具有n 个未 知量的r l 阶线性方程组从中可以求出未知的节点 位移和节点力。 按照上述方法，我们把整个边界离散，将边界划分为 具有n p 个节点的h e 个单元，并在每一个单元上建立局 部坐标系，如图2 5 所示 单元内任意一点的整体坐标、位移和面力分别表示为： r l 7 = 嬉，7 拉j 2 了 l 挑= m ( 善，功”j ，- i 圈2 5 7 。、 l 4 第2 i 页 f 协篮m 。 = 太原理j 二人学碘士学位论文 式中一，“：，f = 1 , 2 ，a ，) 表示单元上第l 节点的坐标、位移 和面力分量，m g ，叩) 为插值函数。边界方程离散化时，单 元类型的选择与划分是影响离散模型误差的重要因素。因 此，对于不同的设计，应选择与之相适应的单元形式来逼 近真实结构。鉴于三维弹性体结构体中，大部分的曲面是 二次曲面，而且处于复杂的应力状态，本文选择了四边形 八节点等参元其插值函数为下式： ( 善，口) 2 j ( 1 + 删+ r ) ( 8 + r 一1 ) ：( 善，r t ) = 一4 - ( 1 一f ) ( 1 + ，7 ) ( 毒一叩+ 1 ) 3 ( 善，7 ) - t 4 0 善) ( 1 一叩埘+ ，7 + 1 ) 神= 丢( 1 + 烈l 一叩) ( 古一，7 - 1 ) 也( 善，7 ) = ( 1 一毒2 ) ( 1 + ，7 ) j 。( 善，叩) ：1 0 一f ) ( 1 - r 2 ) n 7 ( 毒，r t ) = 2 ( 1 一毒2 ) ( 1 一玎) s ( f ，7 ) = 百1 ( 1 + 善) ( 1 - r f l ) 根据上面选择的插值函数，i l i , t 在每个单元内按双- 次曲线变化，当第5 、6 、7 、8 节点有缺时只要将相应 的插值函数取为0 即可即是变节点元的插值函数。 由于离散时要用到单元变换系数j ，其计算公式推导 如下： 单元内任一点的法向矢量；为 太跟理大学坝士学位论文 一 o ro r ，、 舾瓦。百：2 t g - g ：扔】 鲁= c 嘉，詈，嚣， 亳= c 鲁，詈，象， 出，衙1玉，m 1 函3 磊蔷一萌磊 g z = 鲁卺一鲁嚣 岛= 硒a x i 砺a x 2 一善景 则单元变换系数j 为： l ，2 = ( g ? + g ；+ g ；) 依次将每个节点作为单位集中力的作用点即源点，记 作q “( n = l ，2 ，n p ) 。在每个源点上依次旌加单位集中力 毛= l ( j = l ，2 ，3 ) ，并对边界积分方程( 2 17 ) 进行离散化处理 边界积分项化为对每个边界单元的积分之和，于是有： 对物体a ： 一，p ：触。订+ g ：舢。d f + 宁：q 打= t 。封：钮订+ v 二以打+ i ( 一一嚷+ v ) a q ，订 + 掣( 2 1 8 a ) 对物体b 第2 3 页 太腺理1 - 大学蝴十掌位论文 f ，p ：觚。盯+ l 卜吒0 一q a v ；+ 叮：( v ? ，o , n ) 弦+ 眦a v ，+ g + q ：3 ( a v ；一v y 肛= 。卸。a f + l ，( 却? + 吐幻：一屹鲫? ) 订+ f ( 芦v 1 ：v ：_ v ：) 厶叮，d f + 叫( 2 1 8 b ) 上式也可用矩阵形式表示如下： l h “。a “研由5 旬。0 一o “00 一( 蝴曰胁自。+ 自。) i 0 a “00 a h ”一6 “a ；。a ；一仕地8 7 u 茹连一授、 a u 。一 a v “ v ? 5 a v “ 2 a v 。5 ， a u 8 1 。 q 。4 v 厨 i a v “ 2 a q ： l 卸“ - g , - coo o i o = ll l 0 0 。母惨 i3 f a b 。1 + ( 2 1 9 ) 【a 古8j 上式即为第n 步载荷增量迭代的接触方程。其中 a b 。 和 衄8 分别表示a 和b 物体体积力增量项列阵； 疗删= 卜卯4 ，一疗尹，舯m0 “= 【0 p ，一0 尹，铲】。通过联立 迭代求解，可求出边界上的未知的位移和面力。如果采用 全量迭代求解接触问题，则用“，p ，p ，q 分别相应代替 第2 4 页 太原理1 ：大学硕士学位论文 厶”，v ，如，g ，用醒代替讨即可。 对方程( 2 1 9 ) 右端做矩阵乘法方程( 2 1 9 ) 可写成： 【a 】 x = 丑 ( 2 2 0 ) 这就是求解本问题的线性方程组。 2 4 2 奇异积分的处理 在边界积分方程形成的过程中，需要做大量的积分。 而这些积分的精度直接影响线性方程组系数矩阵【a 的精 度，进而影响到位移和应力计算结果的精度。形函数 _ ( 毒。，毒：) 为确定的二次函数，所咀影响积分精度主要是核 函数7 1 ，u 和雅可比值j 在单元中的变化情况。考虑后一情 ，口 况要求我们在边界离散化的过程中注意单元的形态。当单 元为矩形时在单元内i ，( 专，喜，) 为常数，但在实际结构中， 由于形态复杂多变往往很难保证完全正规的单元，这就 要求我们在网格生成的过程中注意单元的形态，使雅可比 值j 在单元中的变化比较小具体的说，边界单元要求不 能有太大的长宽比，棱边的夹角最好为9 0 度，棱边的曲 率不能过大，棱边中节点的位置应尽可能取其中点。还有 一个重要的因素，即核函数在单元中的变化情况，由于u 的主项是1 ， t 。的主项是l r 2 ，为使积分有较好的精度， 要求单元内各点离奇异点的距离相差不能很大。为此，我 们在边界离散的过程中，单元疏密过度要缓慢，边界单元 的长度比不能太大。计算实践说明单元的形态直接影响到 第2 5 碰 太原理t 大学硕士学位论文 结果的精度。 对于上述积分，一般情况下可以采用高斯求积方法， 但当奇异点出现在单元上时，此时被积函数出现奇异性， 直接采用高斯积分方法已经不行。需要进行特殊处理。 对于具有l ，奇异性的积分，我们采用由奇点引出分 割线，将单元划分为二个或三个三角形，在每个三角形中 定义新的局部坐标，这样在新的局部坐标系中的积分，由 于被积函数中出现，。= a ( 亭。，善、) 烈仉，7 7 ，) = o ( r ) ，而使被积函 数的奇异性消除。不失一一般性，我们以奇异点出现在单元 “l i 号和“5 ”号节点为例来说明此问题。 奇异点在“i ”号节点时，如图2 6 所示，边界单元可 以划分为两个三角形，且分别建立新的局部坐标( r ，r ，) 。 4 8 图2 6 奇异点在1 节点图2 7奇异点在5 节点 对三角形1 3 - 4 有： f ：! 二业：型一1 。i 毒2 = ，7 2 = 瓢r = n 0 ) 2o 1 v 2 l = ( 1 u ) 2 川) 1 a ( 印。，) ii 、 ，2 第2 6 页 太琢理j ：大学硕士学位论文 对三角形1 2 - 3 有： 毒i 5 r t 2 f ：! ! 二业：型一i 。2 j i = ( 1 + ，7 2 ) 2 = o ( ，) 奇异点在“5 ”号节点时如图27 所示，边界单元可 以划分为三个三角形： 对三角形5 - 4 1 有： 专2 一( 1 + r ：) 2 芒：! ! 垫兰! 墨! 一l d l2 ( 1 + r 2 ) 4 = o ( r ) 对三角形5 2 3 有： 亭= ( 1 + ，7 ：) 2 f ：翌二丛墨 。2 d 。5 ( 1 + 玎：) 4 = 咐) 对三角形5 3 4 有： 善= 叩。( 1 + 柙2 ) 2 乞2 也 j l = ( 1 + r 2 ) 2 = o ( r ) 同理，对奇异点出现在其他节点情况，我们也不难写 出其坐标变换关系。 对于 阶奇异积分，采用将一组单位刚体位移作为 特解代入边界积分方程中而消除。 第2 7 页 太原理 ：人学硕十学位论文 由于对积分奇异性作了深入的研究和妥善的处理，不 仅产生了计算日和g 矩阵全部元素的方法而且使得数 值计算精度都能得到同样的保证，都可达到高斯积分所具 有的精度。 2 4 3 接区域和类型的判定 在阿文曾经提到，在一定的外载荷作用卜i ，接触边界 上的每一接触点对都必须处于三种接触状态之一，而且满 足接触问题的等式和小等式条件。在接触问题中。由于接 触区域事先无法确定，所以在前处理时，必须假定可能的 接触状态，按接触等式条件建立积分方程计算之后，再 按不等式条件判定，进行迭代求解，接触区域和类型的判 定步骤和准则如下： 一、对本次迭代计算前处于分离状态的点对 1 如果 v ，0 月= = i v ；i v ? n l v ? 1 n l o 则该节点对仍为分离状态。否则该节点对进入接触状 态，必须进行下一步判断。 2 如果 厮百丽万丽占 则该点对进入粘着接触状态( f 为所给定的精度值) ； 否则，进入滑动接触状态。 二、对于本次迭代前处于粘着接触状态的点对 1 如果 一0 笫2 8 负 太原理【大学硕士学位论文 则该节点进入分离状态；否则该节点仍为接触状态， 需进行下一步的判断； 2 如果 矗丽 1 4 研 则进入滑移接触状态：否则，仍为粘着接触状态。 三、对于本次迭代l ； 处于滑移接触状态的点对 1 如果 “0 则该节点进入分离状态；否则，该节点仍为接触状 态，需进行下一步的判断： 2 如果 而了万丽五万s 则该点对进入粘着接触状态( e 为所给定的精度值) ；否 则仍为滑移接触状态。 判定结束后，须根据新的接触状态，对接触区域的有 关数据进行处理，然后进入下一次迭代直到找到证确的 接触区域和接触状态，求得合理的结果。 2 4 5 边界点应力计算 边界上任意点的位移和应力可以由边界元上已知的 节点位移和面力，利用插值关系确定。 设单元e 上任意一点( 毒，7 ) 的位移为“( 善，7 ) ，掌，7 是单 第2 9 页 太原理工大学硕士学位论文 元的局部坐标为简便起见，令f = 善i r l = 乞，则“( 毒，毒：) 对 局部坐标的偏导数为： 加加良 良 吾= 詈= i , i 若 ( i , j = l ，2 ，3k = l ，2 ) 鸭盘，嗨”啦 、 一 现 盘 式中菇和若可以通过插值公式2 - 2 0 的前两式以单元 。t。 e 的节点位移村和节点坐标一表示，均为已知。 联立物理关系和边界平衡条件，得到下列求解边界点 的应力方程： 仃f 一知一 u 枷) 2 0 仃口h j 一， i , j 。i ，2 ，3k 2 1 ，2 出 加 l “= 上 a 。“始t 式中巳是边界点( 毒，毒，) 的应力，盯。和九个位移梯度分量 “共有1 5 个未知量。，由插值公式以节点面力表示。上 j 式中共有1 5 个代数方程，正好用来联立求解得到边界点 应力。 第3 0 页 太原理工大学硕士学位论文 第三章机械构件的形状优化方法 l3 一l 向曩的提出 机械设计中构件的形状优化是比结构尺寸设计更高 层次的设计问愿。用有限元法或边界元法作为分析工具对 构件进行形状优化，使满足工作时对构件功能和性能的要 求，则是近年来发展起来的前沿课题，也是结构优化领域 中目i i 研究得比较少的课题。形状优化与结构尺寸优化的 主要区别在于：刚度矩阵与设计变量之间不一定成i e 比关 系。给不出关于设计变量的解析显式：应力矩阵不是常数 矩阵，而是设计变量的函数。可见，形状优化要比结构尺 寸优化的计算量大得多也困难得多。 尽管有限元法是求解弹性力学问题非常有效的数值 计算方法，但对于形状复杂的构件，比如应力集中问题， 特别是形状剧烈变化处的应力集中问题，为了反映局部区 域的应力变化规律，在该区域必须划分出密集的网格，这 就要求计算机有很大容量。另外，最大应力点一般总是位 于所研究的物体的域边界上，所以外推时又要引起一定的 误差。 由于应力集中多发生于边界，因此就处理弹性力学 的应力集中问题而言，边界元法应当说是一种更有效的数 值计算方法。这种方法只需对边界进行离散，而不必对域 内进行离散，所以给单元离散带来很大方便，它可以直接 得出在域边界上的应力值而无需涉及域内，从而对于计 算应力集中的形状优化问题，由于其有效最大应力点一般 位于域边界上，所以无需计算内部节点的应力，同时，由 第3 i 甄 太原理工大学硬士学位论文 于边界单元的划分简单，所以有利于简化形状优化算法中 的前处理和迭代过程中的单元重划分，因此，用边界元法 更加适合形状优化的特点。 3 2 教学校童 一个具有可变边界l 的弹性连续体，当载荷和位移边 界条件给定时，选择某一形状的边界r 。，使弹性连续体内 应力集中系数获得最小，这样的问题称为求应力集中系数 最小的连续体形状优化问题。r 就是所求的最优解。 求解应力集中系数最小的形状优化问题，目标函数表 达式为： ，2 嘧搿一驯川 ( 3 i ) 式中： f 一应力集中系数。 b 一具有可变边界l 的弹性连续体集合。 r 一具有某一可变边界re r ，的弹性连续体。 q 一弹性体r 内任意点的坐标。 盯( q ) 一点q 处的有效应力值。 p 一定义应力集中系数的标准值。 由于形状改变引起的应力集中最大的点通常是在弹 性体的边界上。因此，在数值计算中，用离散化的边界元 模型来计算连续体表面的应力分布和应力集中系数是很 有利的。本文以指定边界区域上离散节点的有效应力作为 优化应力控制点，以求得该边界上应力集中系数最小的可 变边界形状。一般来讲，指定边界上应力最大的点就是应 第3 2 贞 太原理一夫学硕士学位论文 皇葺j i _ i i 力集中系数的最大点。因此，撅小化指定边界上最大应力 集中系数的形状优化问题，可转化为极小化边界节点上最 大应力的形状优化问题。把目标函数离散化。于是式( 3 1 ) 转化为： f = m n m 似p ( q ) 1 ) l e x 日 式中： ( 3 2 ) x 一设计变量工= ( x ix ：，x ，工。) 1 描述可变边界1 7 ，的 几个参数。 x 一设计变量x 的所有可能取值。 i 一弹性体边界上第1 个节点。 i 一弹性体指定边界上节点的集合。 a ( q ) 一节点i 处有效应力 f 一有效应力。 由于在优化过程中弹性体的边界条件保持不变，所 以各应力控制点的有效应力仅是x 的函数。 进一步，式( 3 2 ) 的极小一极大问题可转化为求x 、 b 厂2 盯( q i ) 一o ，21 ，2 ，3 ， n 。 ( 3 3 ) 式中n 。为边界上应力控制点总数。 在实际形状优化闷题中，由于工程设计需要，同时也 为了避免产生不合理的边界形状，对可变边界需要施加相 应的几何约束，即： 太原理工大学硕士学位论文 g ( 工) so ，朋2 1 , 2 ，3 ，m l 科s i 4 式中： m 一几何约束条件总数。 z ，x 一分别为膏的下界和上界。 对于一般的结构优化问题来说，所用的有效应力是结 构优化中常用的v o nm i s e s 应力： 盯= r 盯2 + 盯2 + 盯2 一仃仃一口盯一盯盯 + 3 c r 2 + 3 0 - ，2