 Advances in Applied Mathematics 应用数学进展, 2013, 2, 74-82 http://dx.doi.org/10.12677/aam.2013.22010 Published Online May 2013 (http://www.hanspub.org/journal/aam.html) Theoretical Inference of the Dynamic Stochastic Finite Element Method Based on the Virtual Displacement Principle Zonghe Zhou, Guangping Yu Naval Deputy Office of 438 Factory, Wuhan Email: zhouzonghe1982@126.com Received Mar. 19th, 2013; revised Mar. 26th, 2013; accepted Apr. 10th, 2013 Copyright © 2013 Zonghe Zhou, Guangping Yu. This is an open access article distributed under the Creative Commons Attribution Li- cense, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract: Virtual displacement principle which is the weak solution integral form of partial differential equa- tions was researched in detail, and axisymmetric isoparametric quadrilateral element with eight-nodes was constructed. Then, based on the virtual displacement principle, the static equilibrium equation of the axisym- metric dynamic problems and the equivalent integral weak form under the boundary conditions of the force were derived, combined with perturbation technology, the dynamic stochastic finite element equation of the axisymmetric isoparametric quadrilateral element with eight-nodes was derived, and the dynamic stochastic finite element method based on the principle of the virtual displacements principle was put forward, the method provides a new solution for researching dynamic response characteristics of engineering objects. Keywords: Virtual Displacement Principle; Perturbation Technology; Axisymmetric Isoparametric Quadrilateral Element with Eight-Nodes; Dynamic Stochastic Finite Element 基于虚位移原理的动态随机有限元方法理论推导 周宗和，余广平 海军驻武汉 438 厂军事代表室，武汉 Email: zhouzonghe1982@126.com 收稿日期：2013 年3月19 日；修回日期：2013 年3月26 日；录用日期：2013 年4月10 日 摘 要：对偏微分方程弱解的积分形式——虚位移原理进行较详细的研究，并构造出八节点四边形轴 对称等参单元。然后从虚位移原理出发，给出轴对称动态问题的平衡方程以及力的边界条件下的等效 积分“弱”形式。结合摄动技术，推导出八节点四边形轴对称等参单元的动态随机有限元方程，提出 了基于虚位移原理的动态随机有限元方法，该方法为研究工程对象的动态响应特性提供了一种新的解 决途径。 关键词：虚位移原理；摄动技术；八节点四边形对称等参单元；动态随机有限元 1. 引言 目前，随机有限元方法[1-7]有多种，这些方法都是以数学、力学分析作为工具，找出结构系统的响应与输入 信号之间的关系，并据此得到结构应力或位移的统计规律，其中包括摄动随机有限元法、Neumann 随机有限元 法和小波随机有限元等。 Copyright © 2013 Hanspub 74  周宗和，余广平 基于虚位移原理的动态随机有限元方法理论推导 首先，工程或物理学中的许多问题，通常是以未知场函数应满足的偏微分方程和边界条件的形式提出来的。 对于工程或物理学中遇到的偏微分方程一般是没有理论解的，即未知函数没有解析表达式。而工程上又需要了 解这些未知函数，所以一般用数值的方法来求解这些偏微分方程。有限元方法就是一种数值求解偏微分方程的 方法，它实际上求解的是偏微分方程的弱解积分形式，所以需要先将偏微分方程变成其弱解积分形式，才能使 用有限元方法。 鉴于此，本文从偏微分方程的弱解积分形式入手，对随机有限元方法进行创新性研究。首先，对偏微分方 程弱解的积分形式——虚位移原理[8]进行较详细的研究，并构造出八节点四边形轴对称等参单元，从虚位移原 理出发，给出轴对称动态问题的平衡方程以及力的边界条件下的等效积分“弱”形式，结合摄动技术，推导出 八节点四边形轴对称等参单元的随机动态有限元方程。 2. 偏微分方程的弱解形式 2.1. 问题的提出 一般地，未知函数 应满足偏微分方程组，表示如下： u 1 20 Au AuA u 在内 (2.1) 域可以是体积域、面积域等，如图 1所示。 0 uA 0 uB 域 x y 0 Figure 1. Solution area 图1. 求解区域 同时未知函数 还应满足边界条件 u 1 20 Bu BuBu 在上 (2.2) 是 域的边界。 要求解的未知函数 u可以是标量函数场(如温度)，也可以是几个变量组成的向量函数场(如位移、应变、应 力等)。, A B是表示对于独立变量(例如空间坐标、时间坐标等)的微分算子。偏微分方程数应和未知场函数的数 目相对应，因此，上述偏微分方程可以是单个的方程，也可以是一组方程。所以在式(2.1)和(2.2)中采用了矩阵 形式。 2.2. 偏微分方程弱解的积分形式——虚位移原理 由于偏微分方程组(2.1)在域 中每一点为零，因此就有 Copyright © 2013 Hanspub 75  周宗和，余广平 基于虚位移原理的动态随机有限元方法理论推导 T 112 2 dVAuAu Au d0 (2.3) 其中 1 2 V (2.4) V是向量函数，称为试验函数或虚位移函数，它是一组和偏微分方程个数相等的任意函数。 式(2.3)是偏微分方程组(2.1)的积分形式。假如 A u是一光滑函数，则可以断言，若积分方程(2.3)对于任意 的V都能成立，则偏微分方程(2.1)必然在域内任一点都得到满足。这个结论的证明是显然的，假如 A u在域内 某些点或—部分子域中不满足，即出现 0Au ，马上可以找到适当的函数V使(2.3)的积分形式亦不等于零。 可见当 A u是一光滑函数时，式(2.3)和(2.1)是等价的。 在很多情况下可以对(2.3)式进行分部积分得到另一种形式 TT dCDu EFu d0 (2.5) 其中 是微分算子，它们中所包含的未知函数导数的阶数比(2.3)式的微分算子,,,CDEF A 低，这样对函数 只需 要求较低阶的连续性就可以了。在(2.5) 式中降低 的连续性要求是以提高 u u 的连续性要求为代价的，由于原来 对 并无连续性要求，但是适当提高对其连续性的要求并不困难，因为它们是可以选择的已知函数。这种降低 对函数 连续性要求的作法在近似计算中，尤其是在有限元方法中是十分重要的。式(2.5)被称为偏微分方程组 (2.1)的弱解积分形式或“弱”形式，或被称为虚位移原理。在使用有限元方法求解前，还需要将(2.2)中的边界 条件代入(2.5)中。 u 3. 构造八节点四边形轴对称等参单元 八节点四边形轴对称等参单元如图2所示，该单元为横截面为八节点四边形的360˚环形等参单元。其横截 面上八个节点的位置坐标 ，各个节点的位移为 ,,1,2,,8 ii rz i ,,1,2,, ii uw i8。 r z Figure 2. Axisymmetric isoparametric quadrilateral element with eight nodes 图2. 八节点四边形轴对称等参单元 如图 2所示，该单元为绕轴的环状单元，在 orz 平面内，八节点四边形轴对称等参单元共有 16 个自由度。 将所有节点上的位移组成一个列阵，记作 ，同样，将所有节点上的各个力也组成一个列阵，记作 ，那 么 z e a e p T 11228 8 e auwuw uw (3.1) T 11228 8 e rzr zrz ppppppp (3.2) 采用的等参变换位移模式和坐标变换式如下： Copyright © 2013 Hanspub 76  周宗和，余广平 基于虚位移原理的动态随机有限元方法理论推导 1) 位移模式 88 11 , iii i ii uNuwN w (3.3) 写成矩阵的形式，得 1 1 2 12 8 2 12 8 8 8 00 0 00 0 e u w u NN N uwNa NN N w u w u (3.4) 2) 应变矩阵 由轴对称问题的几何方程可以推出相应的几何函数矩阵，得 0 10 0 r e z rz r u ruNaB w z zr e a (3.5) 4. 线弹性小变形动态轴对称等参问题的“弱”形式 对于线弹性小变形动态轴对称等参问题，其基本方程[9,10]如下。 如图 2所示，弹性体内任一点沿 方向的平衡方程为 ,rz r rrz r f uu rz r (4.1) rzz rzz f ww rzr (4.2) 其中 ,,, rz rz 为表示体内任一点的应力状态的 4个应力分量， r 为径向正应力， 为环向正应力， z 为轴 向正应力， rz 为径向和轴向之间的剪应力。 , rz f f为单位体积的体积力在 方向的分量，,rz 为密度， 为密度 乘以阻尼系数 。 由平衡方程，得 d d r rrzrz zrz z A A rfu rz rrzr ruuu wwwA fwA (4.3) 其中， ,uw 为 方向上的虚位移。 ,rz 利用分部积分公式将其化为“弱”形式，利用几何方程可得 dd dd r rz zrz rz AA rz rz Al rAruuu rf uf w ArTuTw l wwwA (4.4) Copyright © 2013 Hanspub 77  周宗和，余广平 基于虚位移原理的动态随机有限元方法理论推导 其中， 为单元体的边界。 l 上式表示为向量形式如下 TTTT dddd AA AAl rAru uAru uArufAruTl T d (4.5) 其中 f 为体积力向量 r z f f f (4.6) u 为加速度(位移的二次导数)分量 u uw (4.7) u 为速度(位移的一次导数)分量 u uw (4.8) 将本构方程代入上式，可得 TT TT dd dd AA Al rDAruuu rufAru Tl A (4.9) 将力和位移边界条件代入上式(对于位移边界条件，虚位移,uw 为0)，可得 TT TT dd d AA f Al rDAruuu rufdAru Tl A (4.10) 由此，可得出轴对称等参问题的“弱”形式如下： TT TT dd dd S SS Sf SL rDJSruuuJS rufJSru TL S (4.11) 其中 22 , SL zr JJJ S为坐标变换后求解域的面元，为坐标变换后的单元体的边界。 L 将(3.4)和(3.5) 式代入上(4.11)式，整理后得 TT TT TT TT dd dd eee SS SS ee SfL SL arBDBaJSarNu uJS arNfJSarNTJL (4.12) 5. 轴对称等参问题的动力学随机有限元方程 在许多实际问题当中，影响结构动力学特性的因素，如材料特性，几何特性和载荷特性等都有一定程度随 机性，因此结构的响应 u和 也具有随机性。 为了简化书写，假设任意一个给定的空间连续函数 0 fr ， 0 r为空间的坐标向量，基本随机向量 b的 Copyright © 2013 Hanspub 78  周宗和，余广平 基于虚位移原理的动态随机有限元方法理论推导 联合概率密度函数为 g b， p 为一个小参数。 0 00 00 2 0 d , , Cov,d d jjk jjk pp j kjkjk bEb bgbb fr frb ff ff bbb dbbb b bbb bgbbbb jk 基于摄动理论，将影响结构的随机参数在基本随机变量 b的均值 0 b处作小参数摄动展开，如下： 02 111 02 111 02 111 02 111 0 1 2 1 2 1 2 1 2 nnn pjp jk jjk jjk nnn pjp jk jjk jjk nnn pjp jk jjk jjk nnn pjp jk jjk jjk pj DD DbDbb NN NbNbb BB BbBbb aa ababb ff f 2 111 02 111 02 111 02 111 0 1 1 2 1 2 1 2 1 2 nnn jp jk jk jjk nnn f fpfj pfjk jjk jjk nnn SSpSjj pSjkjk jjk nnn LLp LjjpLjkjk jjk n pjj j bfbb TT TbTbb JJ JbJbb JJ JbJbb rb rr 2 11 02 111 02 111 02 111 02 1 1 2 1 2 1 2 1 2 1 2 nn pjkjk jk nnn pjp jk jjk jjk nnn pjp jk jjk jjk nnn pjj pjkjk jjk n pjjp j rbb uu ububb uu ububb bbb b 11 nn jkj k jk bb (5.1) 将(5.1)式代入(4.12)式，并比较小参数 p 的同次幂系数，可以得到结构摄动格式的各阶随机积分“弱”形式， 其中零阶和一阶随机积分“弱”形式如下： 零阶随机积分“弱”形式： TT TT 00000 000000 TT TT 0000 000 dd dd e ee SS SS ee SfL SL D BaJSarNuuJS fJSarN TJL 0 ar ar 0 0 B N (5.2) Copyright © 2013 Hanspub 79  周宗和，余广平 基于虚位移原理的动态随机有限元方法理论推导 一阶随机积分“弱”形式： T T T T000000 000000000 TTT 00000000 0000000 TT T00000 d ee e jS S jj S e e e SS j j e jS j arBDBaJrB DBaJrBDBaJ rB DBaJrB DBaJrB DBaJS arN uJrN T 00 0000 0 TT 00 000000 TT T T0000 0000 0 000 TT 00 0000 00 TT T0000 d d SjS S SSj j e jS SjS j S SSj j e jS j uJ rNuJ rNuJ rNuJS a rNuJrNuJrNuJ rN uJrN uJS arN fJrNf TT 00 000 000 TTT T T00000000000 0 d d SSSj j S e jfLfLfLfLj j Lj JrN fJrN fJS arN TJrNTJrN TJrN TJL e S Sj (5.3) 假设 和 也有与位移 一样的插值函数，即 u u u e uNa (5.4) e uNa (5.5) 将(5.1)式代入上两式，有 000 000 e e e j j j e j j j uNa uNa uNa uNa (5.6) 由于 是任意的，将(5.6)式代入(5.2)式，经整理，得到单元的零阶动态随机有限元方程： e a 0000000 0 ee ee ee e F e M aCaKaFT (5.7) 其中，为单元的质量矩阵均值，为单元的阻尼矩阵均值， 0e M 0e N 0e K 为单元的线弹性刚度矩阵均值， 为单元的体积力均值， 为单元的面力向量均值，其相应的单元矩阵如下： 0e F 0e F T T 000000 T 000000 T 00000 T 00000 T 00000 2πd 2πd 2πd 2πd 2πd e S S e S S e S S e S S e Ff L 0 L M rN NJS CrNNJS K rBD BJS FrNfJS TrNTJL (5.8) 求解区域离散为 个单元，经组集后可以得到结构的零阶动态随机有限元方程： n Copyright © 2013 Hanspub 80  周宗和，余广平 基于虚位移原理的动态随机有限元方法理论推导 0000000 0 F M aCaKaFT (5.9) 其中，为结构的质量矩阵均值， 为结构的阻尼矩阵均值， 0 M 0 N 0 K 为结构的线弹性刚度矩阵均值， 0 F 为结构的体积力均值， 0 F T为结构的边界面力向量均值。 同理，由(5.3)式，可得到单元的一阶动态随机有限元方程： 00 0 00 ee e ee e jj j ee eee ee F jj jjj Ma CaKa 0 e F TKaMaCa (5.10) 式中， e j M为单元的质量矩阵的变异矩阵， e j C 为单元的阻尼矩阵的变异矩阵， e j K 为单元的线弹性刚度矩 阵的变异矩阵， e j F 为单元的体积力向量的变异向量， e F j T 为单元的边界面力向量的变异向量，分别如下表 示： TT T 000000 000 000 TT 00 000000 2π d e j SSj jj S SSj j S M rNNJ rNNJ rNNJ rNNJ rNNJS (5.11) TT T 00 0000000 000 TT 00 0000 00 2π d e j SSj jj S SSj j CrNNJrN NJrNNJ rNNJ rNNJS S (5.12) T T T 0000 0000 0000 TT 00 0000000 2π d e j SS jj S e SSj j S j K rBDBJrBDBJ rBDBJ rB DBJrB DBaJS T j (5.13) TT T 0000 00000000 2πd ee jS SSS j j j Se F rNfJ rNfJ rNfJ rNfJS (5.14) TT T 0 0000000000 0 2πd ee FjfLfLfL fLj jj Le j TrN TJrNTJrN TJrN TJL T 0 (5.15) 求解区域离散为 个单元，经组集后可以得到结构的一阶动态随机有限元方程： n 00 000 F jj jjj jjj M aCaKaFT KaMaCa (5.16) 6. 小结 本文从另外一个角度对随机有限元方法进行研究，首先对偏微分方程弱解的积分形式—虚位移原理进行较 详细的研究，并根据对八节点四边形等参单元和轴对称问题原有研究的基础上，构造出了八节点四边形轴对称 等参单元。然后从虚位移原理出发，给出轴对称动态问题的平衡方程以及力的边界条件下的等效积分“弱”形 式。结合摄动技术，推导出八节点四边形轴对称等参单元的动态随机有限元方程，提出了基于虚位移原理的动 态随机有限元方法，该方法为研究工程对象的动态响应特性提供了一种新的解决途径。 参考文献 (References) [1] 安伟光, 朱卫兵等. 随机有限元法在不确定性分析中的应用[J]. 哈尔滨工程大学学报, 2002, 23(1): 132-135. Copyright © 2013 Hanspub 81  周宗和，余广平 基于虚位移原理的动态随机有限元方法理论推导 Copyright © 2013 Hanspub 82 [2] 安立强, 王璋奇. 随机约束汽轮机叶片频率的有限元分析[J]. 中国电机工程学报, 2009, 29(2): 95-100. [3] A. Shaker, W. Abdelrahman, et al. Stochastic finite element analysis of the free vibration of functionally graded material plates. Computational Mechanics, 2008, 41(5): 707-714. [4] M. Kaminski. Perturbation-based stochastic finite element method using polynomial response function for the elastic beams. Mechanics Re- search Communications, 2009, 36: 381-390. [5] 安利强. 汽轮机叶片静动特性的随机有限元方法研究[D]. 华北电力大学, 2005. [6] 周宗和, 杨自春等. 基于摄动响应面法的汽轮机转子随机响应特性及灵敏度分析[J]. 汽轮机技术, 2011, 53(4): 41-44. [7] 周宗和, 杨自春等. 汽轮机转子的多变量随机动态响应及可靠性分析[J]. 汽轮机技术, 2011, 53(5): 35-38. [8] 胡海昌. 弹性力学的变分原理及其应用[M]. 北京: 科学出版社, 1981. [9] 王勖成, 邵敏. 有限单元法基本原理和数值方法[M]. 北京: 清华大学出版社, 1997. [10] 曾攀. 有限元方法基本原理[M]. 北京: 清华大学出版社, 2008. |