本发明涉及一种用于获取工程结构静动态力学特性的数值方法， 属于结构力学技术领域。

现有的用于工程结构静动态力学特性分析的技术主要为有限元方 法，它由于具有很好的边界适应性而得到大量和广泛的应用，但目前大 都采用低阶多项式形状函数，因而计算精度不高。当然可以通过两种 途径来提高现行有限元方法的计算精度，第一是通过密化有限元网格 (即h-version方法)，第二是采用高阶次的多项式形状函数(即：p-version 方法和等级有限元方法hierarchical FEM)，近年来，发展了综合前两种方 法的自适应h-p version方法，这些方法虽然可以有效地提高计算精度， 但却使计算量大大地增加，由于所构造的高阶形状函数仍然是采用多 项式(如Legendre正交多项式函数)，因而构造过程复杂、效率很低， 同时，也存在数值稳定性方面的问题。在结构的动态特性分析过程 中，这些缺点则更加突出。

在结构分析中，对于一些情形简单而几何形状又规整的构件，也 可采用经典力学方法，即通过求解微分方程来得到相应的解析解或精 度很高的近似解，这些简单构件包括：杆的拉伸、杆的扭转、梁的弯 曲、方或圆板的平面变形和弯曲变形等，边界条件包括：简支、固 支，由于几何形状及边界条件的严格限制，经典力学方法在实际工程 结构中的应用极为有限。

本发明之目的在于：

1.提出复合单元方法的概念。它可使传统的有限元方法和经典 力学解析解进行复合，从而继承各自方法的优点，具有传统有限元方法 优良的离散特性，即可对任意几何形状的结构进行离散化，又可继承经 典力学解析解的高精度，因而，既具有很好的适应性，又能达到很高的 精度，从而实现高效率。

2.定义两组自由度座标体系来描述复合单元。第一组为节点自 由度座标体系，它主要在于运用传统有限元方法，目的在于对复杂几何 形状进行离散化，第二组为单元的场自由度座标体系(叫作c-自由度)， 它主要在于运用经典力学的解析解，其目的在于提高计算精度，从而实 现高效率。

3.提出“零边界条件”下解析解概念。即由经典力学在零边界 条件下求得解析解，把它作为在复合单元方法中所使用的场自由度座 标体系的基底函数，它是两组自由度座标体系容合的关键，也是传统 有限元法与经典力学解析解复合的基础，本发明申请中可提供适用于 各种边界性要求(如C0、C1问题)的“零边界条件”下解析解。

4.基于所提出的两组自由度座标体系，定义新的形状函数描述(叫 作复合形状函数)，它由两部分组成，其一为传统有限元方法中的形状 函数，其二为经典力学在“零边界条件”下解析解得到的形状函数，同 样，定义新的复合自由度，它也是由节点自由度和场自由度组成。

5.在建立好复合单元的复合形状函数后，可完全按照传统有限元 方法中的过程推导出各种单元(拉伸杆单元、扭转杆单元、梁单元、 平面三节点单元、平面四节点单元、空间四节点四面体单元、空间八 节点六面体单元、三节点弯曲板单元、四节点弯曲板单元......等)的刚 度矩阵及质量矩阵。

6.复合单元方法的实施与传统有限元方法类似，如果将其c-自由 度忽略或取为零，则复合单元方法将退化为传统的有限元方法，也就是 说，传统有限方法是复合单元方法的一种特殊情形。

7.在复合单元方法中，可以通过两种途径来提高计算精度，其一为 与传统有限元方法类似的密化离散网格，(即h-version方法)，其二为增加 c-自由度(叫作c-version方法)，第二种途径是实现复合单元高精度的主 要手段。

本发明用于获取工程结构静动态力学特性的数值复合单元方法， 包括下述步骤： (1)在对结构进行离散化得到各种单元后，定义两组自由度座标体系 来描述单元的位移场，一组为节点自由度座标体系，另一组为场自由 度座标体系(叫作c-自由度)。 (2)基于节点自由度座标体系中的节点自由度q，构造与传统有限元 方法相同的位移场函数UFEM(ξ)及相应的形状函数N(ξ)：

                  UFEM(ξ)＝N(ξ)q   (3)基于场自由度座标体系中的场自由度c，构造由经典力学理论在 “零边界条件”下得到的解析解作为基底的位移场函数UCT(ξ)及相应 的形状函数φ(ξ)：

              UCT(ξ)＝φ(ξ)c     (4)由步骤(2)及(3)得到的位移场函数复合相加而形成复合单元的复合 位移场函数U(ξ)、复合形状函数S(ξ)及复合自由度δ

        U(ξ)＝UFEM(ξ)＋UCT(ξ)

             ＝N(ξ)q＋φ(ξ)c

             ＝S(ξ)δ 其中        S(ξ)＝[N(ξ)φ(ξ)]，δ＝[qc]T (5)在得到复合单元的复合形状函数S(ξ)后，可计算与应变ε(ξ)对应的 几何矩阵B(ξ)

             ε(ξ)＝B(ξ)δ (6)然后计算复合单元的刚度矩阵ke及质量矩阵me $K^e={\int }_V{B}^T\mathrm{DBdV}$ $m^e={\int }_V\rho S^T\mathrm{SdV}$ 其中D为弹性矩阵，ρ为密度。 (7)对结构各单元的刚度矩阵和质量矩阵进行从局部座标到总体座标的 转换，然后再进行组装以得到整体刚度矩阵和整体质量矩阵： $K=\underset{e=1}{\overset{n}{\Sigma }}{T}^{\mathrm{eT}}{k}^e{T}^e$ $M=\underset{e=1}{\overset{n}{\Sigma }}{T}^{\mathrm{eT}}{m}^e{T}^e$ 其中Te为座标转换矩阵。 (8)有了复合单元的整体刚度矩阵和质量矩阵后，可对工程结构进行静 力学分析

                       Kδ＝P 和动力学分析

                       Kδ＝ω2Mδ 其中P为外载荷向量，ω为结构的共振自然频率。 (9)在复合单元方法中，有两种途径来提高数值分析精度，一种是通过 与传统有限元方法类似的密化离散网格，叫作h-version方法，另一 种是增加复合单元的场自由度，即增加c-自由度，叫作c-version方 法，后一种方法的效率远高于前一种。 (10)传统有限元方法可看作为本发明申请的复合单元方法的一种特殊 情况，即将复合单元方法中的场自由度忽略或取为零，则复合单元方 法退化为传统有限元方法。

本发明的方法可用于各种结构的静动力学分析、稳定性分析以及 相应的CAD、FEM软件中。采用复合单元方法可使数值分析效率大大 提高，只需很少的几何离散单元通过增加c-自由度就能获得很高的计算 精度并达到一种超级收敛。

附图说明： 图1表示本发明的实施流程框图。 图2表示复合单元方法中的轴向单元。 图3表示复合单元方法中的扭转单元。 图4表示复合单元法中的弯曲单元。 图5A表示一由7根轴向杆组成的桁架。 图5B表示一与x轴成γ角的任意杆单元。 图6表示一右端固定的悬臂梁。 图7表示各种方法用于计算左端固支轴向拉伸杆第4阶特征频率的相 对误差的比较。 图8表示各种方法用于计算左端固支悬臂梁第4阶特征频率的相对误 差的比较。

下面结合附图详细介绍本发明的内容。

1.对工程结构(或连续体)进行几何离散，如离散为杆单元、梁单 元、板单元、平面三节点单元、空间八节点单元……等。

2.位移场函数的表达，在复合单元方法中，可定义两组自由度座标 体系来描述单元的位移场函数，即

             U(ξ)＝UFEM(ξ)＋UCT(ξ) 其中U(ξ)为单元的复合位移场函数，UFEM(ξ)为基于节点自由度座标体 系的位移场函数，UCT(ξ)为基于单元场自由度座标体系的位移场函数， ξ为无量纲几何位置座标。如果将UFEM(ξ)取为传统有限元方法中的节 点位移场函数，UCT(ξ)取为由“零边界条件”下解析解得到的位移场函 数，则有

                 UFEM(ξ)＝N(ξ)q    

                 UCT(ξ)＝φ(ξ)c 其中N(ξ)为传统有限元方法中的节点形状函数，q为节点自由度，φ(ξ)为 由经典力学的“零边界下解析解”，c为场自由度(即c-自由度)。所以，复 合单元的位移场函数可表达为

             U(ξ)＝N(ξ)q＋φ(ξ)c

                  ＝S(ξ)δ 其中，S(ξ)＝[N(ξ)φ(ξ)]叫作复合单元的复合形状函数，δ＝[qc]T叫作 复合单元的复合自由度。可以看出S(ξ)和δ都是由两个部分组成， 即：基于节点自由度座标体系的传统有限元部分和基于场自由度座标 体系的经典解析解部分。值得注意的是，这里的经典解析解是基于“零 边界条件”下得到的，复合单元方法中所说的“零边界条件”具体可表达 为： 对于C0连续问题，“零边界条件”为

     φr(ξ)|ξ＝0＝0，φr(ξ)|ξ＝1＝0

                                r＝1，2，3，… 其中φr(ξ)为场函数(在工程结构中，φr(ξ)为位移场函数)。 对于C1连续问题，“零边界条件”为

     φr(ξ)|ξ＝0＝0，    φr(ξ)|ξ＝1＝0 $\frac{d{\phi }_r\left(\xi \right)}{\mathrm{d\xi }}{|}_{\xi =0}=0,\frac{{\mathrm{d\phi }}_r\left(\xi \right)}{\mathrm{d\xi }}{|}_{\xi =1}=0$                         r＝1，2，3，… 对于任意的Ck连续问题，“零边界条件”为

       φr(ξ)|ξ＝0＝0，φr(ξ)|ξ＝1＝0 $\frac{d{\phi }_r\left(\xi \right)}{\mathrm{d\xi }}{|}_{\xi =0}=0,\frac{{\mathrm{d\phi }}_r\left(\xi \right)}{\mathrm{d\xi }}{|}_{\xi =1}=0$ $\frac{d^k{\phi }_r\left(\xi \right)}{{\mathrm{d\xi }}^k}{|}_{\xi =0}=0,\frac{d^k{\phi }_r\left(\xi \right)}{d{\xi }^k}{|}_{\xi =1}=0$                            r＝1，2，3，… 以上的“零边界条件”将是获取经典解析解的关键，它将作为基于场自由 度座标体系的形状函数，在大多数情形下，各种单元的φr(ξ)都可以通过 在“零边界条件”下求解微分方程而得到解析表达式。 3.单元刚度及质量矩阵的求取，这一步骤与传统的有限元方法相 同，由于复合单元方法中的形状函数包括有两个部分，因此所得到的单 元刚度及质量矩阵与传统有限元方法中的单元刚度及质量矩阵不同， 后面将具体给出几种类型的单元刚度及质量矩阵。

4.将各单元的刚度及质量矩阵进行装配可形成总体的刚度及质 量矩阵，然后按传统有限元的步骤进行方程求解则可得到结构的位 移、应变、应力、自然频率、振型……等。

下面进一步详述复合单元方法中各种类型单元的位移函数、形状 函数、刚度矩阵以及质量矩阵的表达。 (1)轴向杆单元

图2为一具有轴向位移的杆单元，这是一个C0连续问题，即只要求位 移函数连续，q1和q2表示节点自由度座标体系中的节点自由度， c1c2......cr表示场自由度座标体系中的场自由度(叫作c-自由度)，l为单 元的长度，ξ＝x/l为局部的无量纲几何座标，基于节点自由度座标体 系，轴向位移函数可表达为

         UFEM(ξ)＝q1(1－ξ)＋q2ξ＝N(ξ)q 其中          N(ξ)＝[(1－ξ)ξ]，qT＝[q1q2] 基于场自由度座标体系的位移函数可表达为

UCT(ξ)＝c1φ1(ξ)＋c2φ2(ξ)＋…＋crφ(ξ)＝φ(ξ)c

φi(ξ)，i＝1，2，......r为基底函数，由经典力学理论就杆的拉压问题在“零边 界条件”下所求解而得到，即

             φ1(ξ)＝sin(πξ)

             φ2(ξ)＝sin(2πξ)

                

             φr(ξ)＝sin(rπξ) 所以

     φ(ξ)＝[φ1(ξ)φ2(ξ)…φr(ξ)]

           ＝[sin(πξ)sin(2πξ)…sin(rπξ)]

         cT＝[c1 c2…cr] 那么，杆单元的复合位移场为

U(ξ)＝UFEM(ξ)＋UCT(ξ)

     ＝q1(1－ξ)＋q2ξ＋c1sin(πξ)＋c2sin(2πξ)＋…＋crsin(rπξ)

     ＝S(ξ)δ 其中：

S(ξ)＝[N(ξ)φ(ξ)]

     ＝[(1－ξ)ξsin(πξ)sin(2πξ)…sin(rπξ)] 叫作杆单元的复合形状函数。

          δ＝[qT cT]

            ＝[q1 q2 c1 c2…cr]T 叫作杆单元的复合自由度。

杆单元的轴向应变为 $ϵ=\frac{\partial U\left(x\right)}{\partial x}\text{=}\frac{1}{l}\frac{\partial U\left(\xi \right)}{\partial \xi }$ $=\frac{1}{l}\frac{\partial S\left(\xi \right)}{\partial \xi }•\delta =B\left(\xi \right)•\delta$ 其中 $B\left(\xi \right)=\frac{1}{l}\frac{\partial S\left(\xi \right)}{\partial \xi }$        $=\frac{1}{l}[-11\pi \cos \left(\mathrm{\pi \xi }\right)2\pi \cos \left(2\mathrm{\pi \xi }\right)\dots \mathrm{r\pi }\cos \left(\mathrm{r\pi \xi }\right)]$ 杆单元的刚度矩阵为 $K_L^e={\int }_V{B}^T\mathrm{EBdV}$                   q1  q2   c1  c2   …   …   cr E是杨氏模量，A为杆的横截面积，V为单元的体积。 杆单元的质量矩阵为 $m_L^e={\int }_V\rho S^T\mathrm{SdV}$                  q1   q2     c1 c2    …              …   cr (2)扭转杆单元

图3为一具有扭转位移的杆单元，这是一个C0连续问题，即只要求位 移函数连续，w1及w2表示节点自由度座标体系中的节点自由度， c1、c2、....cr表示场自由度座标体系中的场自由度(叫作c-自由度)，l 为单元的长度，ξ＝x/l为局部的无量纲几何坐标。基于节点自由度座标 体系的扭转位移函数可表达为

        ΘFEM(ξ)＝w1(1－ξ)＋w2ξ＝N(ξ)q 其中：

        N(ξ)＝[(1－ξ)ξ]，q＝[w1 w2]T

基于场自由度座标体系的位移函数可表达为

    ΘCT(ξ)＝c1φ1(ξ)＋c2φ2(ξ)＋…＋crφr(ξ)＝φ(ξ)c 其中φi(ξ)，i＝1，2，.....r为基底函数，它由经典力学理论就杆的扭转问题 在“零边界条件”下求解而得到，即

             φ1(ξ)＝sin(πξ)

             φ2(ξ)＝sin(2πξ)

                    

             φr(ξ)＝sin(rπξ) 所以

      φ(ξ)＝[sin(πξ)sin(2πξ)…sin(rπξ)]

      c＝[c1 c2…cr]     那么，杆单元的复合扭转位移为

     Θ(x)＝ΘFEM(x)＋ΘCT(x)    

          ＝w1(1－ξ)＋w2ξ＋c1sin(πξ)＋c2sin(2πξ)＋…＋cr sin(rπξ)

          ＝N(ξ)q＋φ(ξ)c

          ＝S(ξ)·δ 其中：

    S(ξ)＝[N(ξ)φ(ξ)]

         ＝[(1－ξ)ξsin(πξ)sin(2πξ)…sin(rπξ)] 叫作扭转单元的复合形状函数。

       δ＝[q c]T

         ＝[w1 w2 c1 c2…cn]T 叫作扭转杆单元的复合自由度。

扭转杆单元的剪切应变为 $\mathrm{r\theta x}=y\frac{\partial \Theta \left(x\right)}{\partial x}$ $=\frac{y}{l}\frac{\partial S\left(\xi \right)}{\partial \xi }•\delta =B\left(\xi \right)$ 其中y为距中性轴的距离，B(ξ)为 $B\left(\xi \right)=\frac{y}{l}\frac{\partial S\left(\xi \right)}{\partial \xi }$ $=\frac{y}{l}[-11\pi \cos \left(\text{πξ}\right)2\pi \cos \left(2\mathrm{\pi \xi }\right)\dots \mathrm{r\pi }\cos \left(\mathrm{r\pi \xi }\right)]$ 扭转杆单元的刚度矩阵为 $K_T^e={\int }_V{B}^T\mathrm{GBdV}$

              w1  w2    c1 c2  …     …   cr G为弹性剪切模量，J为扭转杆横截面的极惯性矩，V为单元的体积。

扭转杆单元的质量矩阵为 $m_T^e={\int }_V\rho S^T\mathrm{SdV}$

                w1   w2      c1 c2    …           …    cr (3)弯曲梁单元

图4为一具有横向位移的弯曲梁单元，这是一个C1连续问题，即要 求挠度的一阶导数连续，v1和v2表示节点自由度座标体系中的节点横 向挠度位移，θ1及θ2表示节点自由度座标体系中的节点弯曲转角， c1，c2，....，cr表示场自由度座标体系中的场自由度(叫做c-自由度)，l为 单元的长度，ξ＝x/l为局部的无量纲几何座标，基于节点自由度座标体 系的弯曲，挠度函数可表达为 WFEM(ξ)＝v1(1－3ξ2＋2ξ3)＋θ1l(ξ－2ξ2＋ξ3)＋v2(3ξ2－2ξ3)＋θ2l(ξ3－ξ2)

    ＝N(ξ)q 其中：N(ξ)＝[(1－3ξ2＋2ξ3)l(ξ－2ξ2＋ξ3)(3ξ2－2ξ3)l(ξ3－ξ2)]

  q＝[v1 θ1 v2 θ2]T 基于场自由度座标体系，弯曲挠度位移函数可表达为

WCT(ξ)＝c1φ1(ξ)＋c2φ2(ξ)＋…＋crφr(ξ)＝φ(ξ)c φi(ξ)，i＝1，2，.....r为基底函数，它由经典力学理论就梁的弯曲问题在“零 边界条件”下求解而得到，即 ${\phi }_1\left(\xi \right)=F_1({\lambda }_1^{*},\xi )$ ${\phi }_2\left(\xi \right)=F_2({\lambda }_2^{*},\xi )$                                      ${\phi }_r\left(\xi \right)=F_r({\lambda }_r^{*},\xi )$ 其中： $F_i({\lambda }_i^{*},\xi )={\sin \lambda }_i^{*}\xi -{\sinh \lambda }_i^{*}\xi -\left(\frac{\sin {\lambda }_i^{*}-\sinh {\lambda }_i^{*}}{\cos {\lambda }_i^{*}-\cosh {\lambda }_i^{*}}\right)(\cos {\lambda }_i^{*}\xi -\cos {\lambda }_i^{*}\xi )$                             i＝1，2，3，… 并且λi*满足关系： ${\cos \lambda }_i^{*}•\cos {\mathrm{h\lambda }}_i^{*}=1$                       i＝1，2，3，… 即： ${\lambda }_1^{*}=40730041$ ${\lambda }_2^{*}=7.853205$ ${\lambda }_3^{*}=10.995608$ ${\lambda }_4^{*}=14.137165$                                        ${\lambda }_r^{*}=(r+0.5)\pi ,$                     r＞4 所以： $\phi \left(\xi \right)=[F_1({\lambda }_1^{*},\xi )F_2({\lambda }_2^{*},\xi )\dots F_r({\lambda }_r^{*},\xi )]$ c=[c1 c2…cr]T 那么，梁单元的复合挠度位移为 W(ξ)＝WFEM(ξ)+WCT(ξ)

 ＝v1(1－3ξ2＋2ξ3)+θ1l(ξ－2ξ2+ξ3)＋v2(3ξ2－2ξ3)＋θ2l(ξ3－ξ2)＋       $c_1{F}_1({\lambda }_1^{*},\xi )+c_2{F}_2({\lambda }_2^{*},\xi )+\dots +c_r{F}_r({\lambda }_r^{*},\xi )$      ＝N(ξ)q＋φ(ξ)c      ＝S(ξ)·δ 其中：

 S(ξ)＝[N(ξ)φ(ξ)]

     ＝[(1－3ξ2＋2ξ3)l(ξ－2ξ2＋ξ3)(3ξ2－2ξ3)l(ξ3－ξ2) $F_1({\lambda }_1^{*},\xi )F_2({\lambda }_2^{*},\xi )\dots F_r({\lambda }_r^{*},\xi )]$ 叫作梁单元的复合形状函数。

    δ＝[qT cT]

      ＝[v1 θ1 v2 θ2 c1 c2…cr]T 叫作梁单元的复合自由度。

梁单元的应变为： $ϵ=-\stackrel{-}{y}\frac{{\partial }^{2}W\left(x\right)}{\partial x^2}$ $=-\frac{\stackrel{-}{y}}{l^2}\frac{{\partial }^{2}S\left(\xi \right)}{{\partial \xi }^2}•\delta =B(\xi •\delta )$ 其中y为距中性轴距离，B(ξ)为 $B\left(\xi \right)=-\stackrel{-}{y}[\frac{1}{l^2}(12\xi -6)\frac{1}{l}(6\xi -4)-\frac{1}{l^2}(12\xi -6)\frac{1}{l}(6\xi -2)$ $F_1\text{'}\text{'}{(\lambda }_1^{*},\xi )F_2\text{'}\text{'}({\lambda }_2^{*},\xi )\dots F_r\text{'}\text{'}({\lambda }_r^{*},\xi )]$ $F_r\text{'}\text{'}({\lambda }_r^{*},\xi )=$ $-{\lambda }_r^{*2}[\sin {\lambda }_r^{*}\xi +\sinh {\lambda }_r^{*}\xi -\left(\frac{{\sin \lambda }_r^{*}-\sinh {\lambda }_r^{*}}{\cos {\lambda }_r^{*}-\cosh {\lambda }_r^{*}}\right)\cos {\lambda }_r^{*}\xi +\cosh {\lambda }_r^{*}\xi$                            r＝1，2，3，… $\xi =\frac{x}{l}$ ${\lambda }_1^{*}=4.730041$ ${\lambda }_2^{*}=7.853205$ ${\lambda }_3^{*}=10.995608$ ${\lambda }_4^{*}=14.137165$ ${\lambda }_r^{*}=(r+0.5)\pi ,$      r＞4 弯曲梁单元的刚度矩阵为 $k_B^e={\int }_V{B}^T\mathrm{DBdV}$ $=E{\int }_A{\int }_0^l{B}^T\mathrm{BdxdA}$

                  v1  θ1   v2 θ2       c 其中，I为梁单元横截面的惯性矩，V为单元体积，子矩阵kBcc为

           c1           c2        c3        c4           c5    …  cr 弯曲梁单元的质量矩阵为 $m_B^e={\int }_V\rho S^T\mathrm{SdV}$ $=\mathrm{\rho A}{\int }_0^l{S}^T\mathrm{Sdx}$

                    v1    θ1    v2    θ2        c ·V为单元体积，子矩阵mBqc为

           v1         θ1          v2         θ2 子矩阵mBcc为

            c1        c2       c3       c4    c5 …  cr (4)其它的复合单元

对于平面三节点单元、平面四节点单元、空间四节点四面体单 元、空间八节点六面体单元、三节点弯曲板单元、四节点弯曲板单元 等各种类型的复合单元以及对应的等参复合单元，完全可以按照前面 所述的方法构造其相应的基于节点自由度座标体系和场自由度座标体 系的复合位移函数、复合形状函数、应变、刚度矩阵以及质量矩阵。 (5)应用实例

图5A表示一个由7根杆组成的桁架。现分析它的各阶共振频率和 模态，图中的参数如下：L＝2l＝4m，h＝2m，横截面积A＝0.001m2，杆的质 量密度ρ＝8000Ns2/m4，杨氏模量E＝2.1*1011N/m2。

图5B为一与x轴成γ角的任意杆单元，[qi qj]为节点的局部轴向位 移，[ui vi uj vj]为整体节点位移，它们之间的转换关系为

         qi＝uicosγ＋visinγ

         qj＝ujcosγ＋vjsinγ 取任意复合杆单元的局部“复合自由度”δe为

        δe＝[qi qj c1 c2…cr]T 取任意复合杆单元的整体“复合自由度” δe为

         δe＝[ui vi uj vj c1 c2…cr]T 那么δe与 δe之间的转换关系为

          δe＝Te δe 其中：

         ui   vi   uj    vj    c1 c2 …    …   cr 那么，在整体座标体系中，杆单元的刚度矩阵为

Ke＝TeTkeTe

       ui            vi        uj       vj    c1…cr 杆单元的质量矩阵为

Me＝TeTmeTe

         ui           vi         uj        vj      c1 …  cr 分别对图5A中桁架中的7个杆单元形成各自的基于整体座标的单元刚 度矩阵和质量矩阵，然后将它们组装，形成总刚度矩阵和质量矩阵 $K=\underset{e=1}{\overset{7}{\Sigma }}{K}^e$ $M=\underset{e=1}{\overset{7}{\Sigma }}{M}^e$ 再求解以下振动方程

         Kδ＝ω2Mδ 如果每个杆单元只取2个场自由度(c-自由度)，则求解的前10阶自然频 率为

ω1＝1648.26Hz，ω2＝1741.32Hz，

ω3＝3113.83Hz，ω4＝4567.69Hz，

ω5＝4829.70Hz，ω6＝7379.96Hz，

ω7＝7532.30Hz，ω8＝8047.93Hz，

ω9＝9997.48Hz，ω10＝10567.43Hz， 如果用传统的有限元方法求解此问题，则得到的前6阶自然频率为

ω1＝1683.52Hz，ω2＝1776.28Hz，

ω3＝3341.37Hz，ω4＝5174.35Hz，

ω5＝5678.1 8Hz，ω6＝8315.40Hz， 下面再给出一个梁单元应用的例子。

图6所示为一右端固定的悬臂梁，L为梁的长度，ρ为密度，E为杨氏 模量，I为横截面的惯性矩，现分析它的各阶共振频率和模态。

如果只用一个复合梁单元来处理这一问题，取4个场自由度(c-自由 度)，它的总刚度矩阵K为

       v1 θ1         c1        c2        c3         c4 总质量矩阵为

           v1      θ1           c1    c2      c3        c4 令 ${\lambda }_i^4=\frac{{\mathrm{\rho L}}^4}{\mathrm{EI}}{{\omega }_i}^2$                         i＝1，2，… 其中ωi为共振自然频率。

由复合单元方法计算出的各阶λi如下：

λ1＝1.875109，    λ2＝4.694419，

λ3＝7.857543，    λ4＝11.00451， 由传统有限元方法取3个单元计算出的各阶λi如下：

λ1＝1.875199，    λ2＝4.701793，

λ3＝7.903542，    λ4＝11.86048， 由解析解得到的各阶精确的λi如下：

λ1＝1.875104，    λ2＝4.694091，

λ3＝7.854757，    λ4＝10.99554，

图7为各种方法用于计算左端固定拉伸杆的第4阶特征频率(λ4)的 误差比较，图8为各种方法用于计算左端固定悬臂梁的第4阶特征频率 (λ4)的误差比较，可以看出，在相同计算量的前提下，复合单元方法 的计算结果的误差大大低于常规有限元方法的误差，其计算效率高出 10倍左右，取较多场自由度(c-自由度)的结果明显优于取较少场自 由度的结果。