6-1 6-1 基本量及基本方程的矩阵表示基本量及基本方程的矩阵表示6-2 6-2 有限单元法的概念有限单元法的概念 6-3 6-3 单元的位移模式与解答的收敛性单元的位移模式与解答的收敛性 6-4 6-4 单元的应变列阵和应力列阵单元的应变列阵和应力列阵6-5 6-5 单元的结点力列阵与劲度列阵单元的结点力列阵与劲度列阵6-6 6-6 荷载向结点移置荷载向结点移置 单元的结点荷载列阵单元的结点荷载列阵6-7 6-7 结构的整体分析结构的整体分析 结点的平衡方程组结点的平衡方程组6-8 6-8 解题的具体步骤解题的具体步骤 单元的划分单元的划分6-9 6-9 计算成果的整理计算成果的整理6-10 6-10 计算实例计算实例第六章第六章 用有限单元法解平面问题用有限单元法解平面问题有限单元法工程应用实例1头盔撞击试验仿真模型与结果头盔撞击试验仿真模型与结果 有限单元法工程应用实例2 高强钢板厚度高强钢板厚度10mm，材料考虑应变率影响和失效，不受任何约束。模，材料考虑应变率影响和失效，不受任何约束。模拟受初始速度为拟受初始速度为120m/s和和180m/s钢球的冲击过程。钢球的冲击过程。 穿甲试验仿真穿甲试验仿真 初速度为120m/s初速度为180m/s有限单元法工程应用实例3 动画显示的是地动画显示的是地基中心点的沉降基中心点的沉降随线性荷载的变随线性荷载的变化过程，云图显化过程，云图显示莫尔库仑材料示莫尔库仑材料的塑性区形成和的塑性区形成和大变形塑性流动大变形塑性流动过程。过程。 群桩复合地基承载力计算结果群桩复合地基承载力计算结果6-16-1 基本量及基本方程的矩阵表示基本量及基本方程的矩阵表示体力列阵：体力列阵：面力列阵：面力列阵：应力列阵：应力列阵：应变列阵：应变列阵：位移列阵：位移列阵： 物理方程：物理方程：称为称为弹性矩阵弹性矩阵。对于对于平面应变问题平面应变问题，只需将，只需将弹性矩阵弹性矩阵 D D 中的中的E E、 分别分别换成换成 即可。即可。 （平面应力问题）（平面应力问题）则则虚功方程虚功方程可用矩阵表示为：可用矩阵表示为：几何方程：几何方程：此外，用限单元法还要用到此外，用限单元法还要用到虚功方程虚功方程：现将现将虚位移虚位移及与该虚位移相应的及与该虚位移相应的虚应变虚应变表示为：表示为：对连续变形体，它可以代替平衡微分方程和应力边界条件。对连续变形体，它可以代替平衡微分方程和应力边界条件。1 1、对连续体进行、对连续体进行离散化离散化。6-2 6-2 有限单元法的概念有限单元法的概念 有限单元法是有限单元法是用由有限多个、有限大小的单元在有限个结点相互连接的用由有限多个、有限大小的单元在有限个结点相互连接的集合体来近似原来的连续体集合体来近似原来的连续体，当上述单元足够小从而划分网格足够密时，就，当上述单元足够小从而划分网格足够密时，就可以真实地模拟原连续体。可以真实地模拟原连续体。有限单元法分析的基本步骤：有限单元法分析的基本步骤：2 2、单元分析单元分析：（1 1）选择适当的）选择适当的位移模式位移模式，用单元结点位移（为基本未知量）来表示单元，用单元结点位移（为基本未知量）来表示单元 内任一点的位移，即要建立如下关系式内任一点的位移，即要建立如下关系式: : 对于平面问题，最简单而对于平面问题，最简单而常用的单元是常用的单元是三角形单元三角形单元。在。在平面应力问题中，它们是三角平面应力问题中，它们是三角板，在平面应变问题中，它们板，在平面应变问题中，它们是三棱柱。是三棱柱。结点结点铰接点铰接点 d d e 称为称为单元结点位移列阵单元结点位移列阵。 （2）应用几何方程，求出单元的应变，即：）应用几何方程，求出单元的应变，即：（3）应用物理方程，求出单元的应力，即：）应用物理方程，求出单元的应力，即：其中其中 S 称为称为应力转换矩阵。应力转换矩阵。 ijmuiviumvmyxujvjO其中其中 N 称为称为形函数矩阵形函数矩阵。其中其中 B 称为称为应变转换矩阵。应变转换矩阵。 （5）将作用在单元上的外荷载按）将作用在单元上的外荷载按虚功相等虚功相等的原则，的原则， 移置到单元各结点处，成为移置到单元各结点处，成为单元结点荷载单元结点荷载：（4）由于单元产生了应力，则在单元的边界及内部作用有与之平衡的面力和）由于单元产生了应力，则在单元的边界及内部作用有与之平衡的面力和 体力；现将其按体力；现将其按虚功相等虚功相等的原则移置到单元各个顶点处，作为结构其它的原则移置到单元各个顶点处，作为结构其它 部分通过结点对此单元的作用力，称部分通过结点对此单元的作用力，称单元结点力单元结点力，再利用，再利用虚功方程虚功方程，得：，得：即为单元结点力，即为单元结点力， k 称为称为单元劲度矩阵单元劲度矩阵。对各结点进行平衡分析，列平衡方程并组集对各结点进行平衡分析，列平衡方程并组集得到得到整体结点平衡方程组整体结点平衡方程组：其中其中 K 称称整体劲度矩阵整体劲度矩阵。ijmOxFixFiyyFmyFmxFjyFjx其中：其中：3 3、整体分析整体分析：6-3 6-3 单元的位移模式与解答的收敛性单元的位移模式与解答的收敛性对对三结点三角形单元三结点三角形单元，假设位移分量只是坐标的线性函数，即：，假设位移分量只是坐标的线性函数，即： 由左边三个方程求解由左边三个方程求解a a1 1 、 a a2 2 、 a a3 3 ，右边三个方程求解，右边三个方程求解 a a4 4 、 a a5 5 、 a a6 6 。再代回。再代回u、v 中，得中，得 在在i、j、m三个结点，位移应当等于结点位移，即：三个结点，位移应当等于结点位移，即：ijmuiviumvmyxujvjO一、位移模式：一、位移模式：其中，其中，Ni 、Nj、 Nm 称称形函数形函数，其表达式为，其表达式为为单元为单元 i j m 的面积。为使面积不致为负，在图示的面积。为使面积不致为负，在图示坐标系中坐标系中i j m 的次序须是的次序须是逆时针逆时针的。的。而而ijmuiviumvmyxujvjO分别为分别为系数行列式系数行列式第一、二、三列各元素的代数余子式第一、二、三列各元素的代数余子式。则单元内任一点的位移可用矩阵表示为：则单元内任一点的位移可用矩阵表示为：为为单元结点位移列阵单元结点位移列阵。为为形函数矩阵形函数矩阵。其中：其中：简写为简写为二、形函数的几何意义及性质：二、形函数的几何意义及性质：记三角形单元记三角形单元 i j m 内的任一点为内的任一点为P (x , y)，则知形函数的几何意义为：，则知形函数的几何意义为：ijmuiviumvmyxujvjO由此几何意义容易看出形函数具有如下性质由此几何意义容易看出形函数具有如下性质： 为了保证有限单元法解答的收敛性，必须使位移模式能够正确反映物体为了保证有限单元法解答的收敛性，必须使位移模式能够正确反映物体的真实位移形态，具体说来，就是要满足下列的真实位移形态，具体说来，就是要满足下列三方面的条件三方面的条件。 位移模式必须能反映单元的位移模式必须能反映单元的刚体位移刚体位移。 位移模式必须能反映单元的位移模式必须能反映单元的常量应变常量应变。 位移模式应当尽可能反映位移模式应当尽可能反映位移的连续性位移的连续性。在在 i j 及及i m 两边的中点，两边的中点，在三角形在三角形 i j m 的形心的形心，ijmuiviumvmyxujvjO三、解答的收敛性：三、解答的收敛性： 解答的收敛性是指：当单元的尺寸逐步取小时有限单元法的解答收敛解答的收敛性是指：当单元的尺寸逐步取小时有限单元法的解答收敛于真实的解答。于真实的解答。注意注意： 为必要条件，为必要条件， 为充分条件。为充分条件。6-4 6-4 单元的应变列阵和应力列阵单元的应变列阵和应力列阵其中其中 B 称称应变转换矩阵应变转换矩阵，可写成，可写成或简写为或简写为将位移将位移u、v 代入几何方程，可得用结点位移表示的单元应变：代入几何方程，可得用结点位移表示的单元应变：再将单元的应变代入物理方程，得到用结点位移表示的单元应力再将单元的应变代入物理方程，得到用结点位移表示的单元应力其中其中：可简写为可简写为：称称应力转换矩阵应力转换矩阵。可写成分块形式可写成分块形式注意：注意：由于矩阵由于矩阵 B 的元素都是常量，可见应变的元素都是常量，可见应变 e e 的元素也是常量。因此的元素也是常量。因此三结点三角形单元也称为三结点三角形单元也称为平面问题的常应变单元平面问题的常应变单元。（平面应力）（平面应力）注意：注意：在每个单元中，应力分量也是常量。由于相邻单元一般将具有不同的在每个单元中，应力分量也是常量。由于相邻单元一般将具有不同的应力，因而在它们的公共边上，应力，因而在它们的公共边上，应力并不连续应力并不连续。6-5 6-5 单元的结点力列阵和劲度矩阵单元的结点力列阵和劲度矩阵另一方面，由另一方面，由虚功方程虚功方程有有ijmOxFixFiyyFmyFmxFjyFjx 由于单元产生了应力，则在单元的内部及边界作用有与之平衡的体力和由于单元产生了应力，则在单元的内部及边界作用有与之平衡的体力和面力；现将其按面力；现将其按虚功相等虚功相等的原则移置到单元各个顶点处，成为的原则移置到单元各个顶点处，成为单元结点力单元结点力：设单元结点设单元结点 i 、j 、m 发生了虚位移，即：发生了虚位移，即：则有：则有：于是由以上两式，得于是由以上两式，得（ t 单元厚度单元厚度 ）则：则：其中其中 k 称为称为单元劲度矩阵单元劲度矩阵。由于虚位移可以是任意的，再令：由于虚位移可以是任意的，再令：代入上式：代入上式：注意到：注意到： 对三结点三角形单元对三结点三角形单元 i j m ， k 可写成分块形式：可写成分块形式：（平面应力）（平面应力）例：图示等腰直角三角形单元例：图示等腰直角三角形单元 i j m 。试写出单元。试写出单元 的应力转换矩阵的应力转换矩阵 S 和劲度矩阵和劲度矩阵 k 。在图示坐标下在图示坐标下解：解：ijmxyaa下面求单元的劲度矩阵下面求单元的劲度矩阵 k ：单元劲度矩阵单元劲度矩阵 k 元素的元素的力学意义力学意义：同理可求其它分块。最后得：同理可求其它分块。最后得：单元劲度矩阵单元劲度矩阵 k 的的特点特点：注意到只有结点注意到只有结点 i 上有位移，其它结点位移均为零，则有上有位移，其它结点位移均为零，则有例：上例中若例：上例中若 j、m 处为固定铰链支座，处为固定铰链支座， 结点结点 i 处作用有处作用有 水平力水平力P 和竖直力和竖直力 P ，求单元位移分量，求单元位移分量 u 、v 。解：知：解：知：其中：其中：由于：由于：即：即：于是于是所以所以而由上例结果知而由上例结果知：解得解得6-6 6-6 荷载向结点移置荷载向结点移置 单元的结点荷载列阵单元的结点荷载列阵 本节将按本节将按虚功相等原则虚功相等原则将单元所受实际荷载（集中荷载、体力、面力）将单元所受实际荷载（集中荷载、体力、面力）向单元结点移置而成为向单元结点移置而成为单元结点荷载单元结点荷载。ijmFLixFLiyFLmyFLmxFLjyFLjxfPyfPxMxOy在一定位移模式下，这样移置的结果是唯一的，且原荷载与移置后的结点在一定位移模式下，这样移置的结果是唯一的，且原荷载与移置后的结点荷载在向同一点简化时，具有相同的主矢量和主矩。荷载在向同一点简化时，具有相同的主矢量和主矩。一、集中荷载一、集中荷载：设在点设在点 M 有集中力有集中力：为单位厚度上力的大小。则可推得为单位厚度上力的大小。则可推得：即即：（ t 为厚度）为厚度）ijmFLixFLiyFLmyFLmxFLjyFLjxfyfxxOy二、分布体力二、分布体力：设单元设单元受有分布体力受有分布体力：则可将微分体积则可将微分体积 t dx dy 上的体力上的体力：即即：当作集中荷载，利用以上结果并积分得到当作集中荷载，利用以上结果并积分得到：例：例：设单元设单元 i j m 的密度为的密度为 r r ，试求自重的等效结点，试求自重的等效结点荷载。荷载。解：解：由于由于则则ijmFLixFLiyFLmyFLmxFLjyFLjxxOy三、分布面力三、分布面力：设单元某边上设单元某边上受有分布面力受有分布面力：则可将微分面积则可将微分面积 t ds 上的面力上的面力：即即：当作集中荷载，利用其结果并积分得到当作集中荷载，利用其结果并积分得到：例：例：设单元在设单元在 i j 边上受有沿边上受有沿 x 方向的均布面力方向的均布面力 q ，试求等效结点，试求等效结点荷载。荷载。结论：结论：在采用在采用线性位移模式线性位移模式的情况下，单元荷载向结点的移置与理论力学的情况下，单元荷载向结点的移置与理论力学的的刚体静力等效原则刚体静力等效原则完全一致。完全一致。解：解：由于由于则则6-7 6-7 结构的整体分析结构的整体分析 结点的平衡方程组结点的平衡方程组 由前面的单元分析知，在有限单元法中，由前面的单元分析知，在有限单元法中，各个单元只受单元结点力各个单元只受单元结点力 F e作用且处于平衡状态作用且处于平衡状态。本节则进一步说明如何进行结点平衡分析。本节则进一步说明如何进行结点平衡分析。 对于结构中的任一结点对于结构中的任一结点n ，作用于其上的力有两种：一种是围绕结点，作用于其上的力有两种：一种是围绕结点 n的单元对结点的单元对结点 n 的作用力的作用力 Fn ，这种力为，这种力为单元结点力的反作用力单元结点力的反作用力；另一种；另一种是作用于结点是作用于结点 n 上的上的整体结点荷载整体结点荷载 FLn ，包括从围绕结点，包括从围绕结点 n的单元上移置的单元上移置过来的单元结点荷载，以及本来就作用于结点过来的单元结点荷载，以及本来就作用于结点 n 上的集中荷载或支反力。上的集中荷载或支反力。则结点则结点n的平衡方程为：的平衡方程为：其中其中 是对围绕结点是对围绕结点n的单元求和。的单元求和。再代入单元结点力与单元结点位移的关系：再代入单元结点力与单元结点位移的关系：并将同一结点位移的项合并，可得：并将同一结点位移的项合并，可得：( ( q 为结点总数为结点总数 ) 称为称为整体结点位移列阵整体结点位移列阵称为称为整体结点荷载列阵整体结点荷载列阵 K 为为整体劲度矩阵整体劲度矩阵，它是由单元劲度矩阵按结点的局部编码（，它是由单元劲度矩阵按结点的局部编码（ i , j , m ）与整体编码的对应关系组集而成的。与整体编码的对应关系组集而成的。结构整体分析的步骤：结构整体分析的步骤：1、组集整体劲度矩阵、组集整体劲度矩阵 K ；3、引入位移约束条件；、引入位移约束条件；4、求解整体结点平衡方程组、求解整体结点平衡方程组 K d d = FL ，得结点位移，得结点位移 d d ；5、求各单元的位移、应力。、求各单元的位移、应力。将所有结点的平衡方程按结点整体编码组集，得将所有结点的平衡方程按结点整体编码组集，得整体结点平衡方程组整体结点平衡方程组：其中：其中：2、组集整体结点荷载列阵、组集整体结点荷载列阵 FL ； 如图正方形薄板划分为两个单元，厚度为如图正方形薄板划分为两个单元，厚度为t t，密度为密度为r r，弹性模量为，弹性模量为E E ，取泊松比，取泊松比 = 0 = 0 ，各，各单元直角边长为单元直角边长为a , ,结点的局部编码与整体编码结点的局部编码与整体编码的对应关系如下：的对应关系如下：单元号单元号局部编码局部编码整体编码整体编码i41j14m32要求各结点的位移和各单元的应力。要求各结点的位移和各单元的应力。下面用实例说明以上求解过程。下面用实例说明以上求解过程。一、组集整体劲度矩阵一、组集整体劲度矩阵 K 整体劲度矩阵整体劲度矩阵 K 以整体编码排列，而单元以整体编码排列，而单元劲度矩阵劲度矩阵 k e以局部编码排列，故要把单元劲度以局部编码排列，故要把单元劲度矩阵的各子矩阵按对应的整体编码集成。矩阵的各子矩阵按对应的整体编码集成。单元号单元号局部编码局部编码整体编码整体编码i41j14m32组集整体劲度矩阵组集整体劲度矩阵 K 的方法为：的方法为： 将将 K 的全部元素充零；的全部元素充零； 逐个单元地建立逐个单元地建立 k e ，然后按单元，然后按单元 e 中中局部编码与整体编码的对应关系，局部编码与整体编码的对应关系， 对所有单元完成上述叠加后，就形成了整体劲度矩阵对所有单元完成上述叠加后，就形成了整体劲度矩阵 K 。整体劲度矩阵整体劲度矩阵 K 的的特点特点：对称对称对角线元素对角线元素0奇异。奇异。二、组集整体结点荷载列阵二、组集整体结点荷载列阵 FL 结点结点 n 上的上的整体结点荷载整体结点荷载 FLn ，包括从围绕结点，包括从围绕结点 n的单元上移置过来的单元上移置过来的单元结点荷载，以及本来就作用于结点的单元结点荷载，以及本来就作用于结点 n 上的集中荷载或支反力。上的集中荷载或支反力。各单元的结点荷载列阵为：各单元的结点荷载列阵为： 本来就作用于各结点上的本来就作用于各结点上的集中荷载或支反力为：集中荷载或支反力为：组集关系为：组集关系为：最后得结构的整体结点荷载列阵为：最后得结构的整体结点荷载列阵为：三、引入位移约束条件三、引入位移约束条件结构的位移约束（即已知的结点位移分量）为：结构的位移约束（即已知的结点位移分量）为：于是，于是，整体结点位移列阵整体结点位移列阵简化为简化为引入位移约束的方法为：引入位移约束的方法为：（1）将与已知的结点位移分量相应的平衡方程）将与已知的结点位移分量相应的平衡方程 去掉，从而去掉了去掉，从而去掉了 FL 中的未知支反力。中的未知支反力。 本例中即划去本例中即划去 K 、 FL 中的第中的第1、2、3、5、6、8 行。行。（2）将零结点位移分量代入其余的平衡方程，）将零结点位移分量代入其余的平衡方程， 从而去掉了从而去掉了 d d 中的已知结点位移分量。中的已知结点位移分量。 本例中即划去本例中即划去 K 中的第中的第1、2、3、5、6、8 列和列和 d d 中的零结点位移分量。中的零结点位移分量。整体结点荷载列阵整体结点荷载列阵简化为：简化为：（为何不可立即求解方程组（为何不可立即求解方程组 K d d = FL ？）？）整体劲度矩阵整体劲度矩阵简化为：简化为：最后得：最后得：可见引入位移约束条件后，可见引入位移约束条件后， K 仍为对称矩阵，但却是非奇异的。仍为对称矩阵，但却是非奇异的。四、求解整体结点平衡方程组四、求解整体结点平衡方程组 K d d = FL ，得结点位移，得结点位移 d d 解得解得 五、求各单元应力五、求各单元应力问题：如何求支座的约束反力？问题：如何求支座的约束反力？(4)列出结构的整体结点平衡方程组列出结构的整体结点平衡方程组 K d d = FL ，引入位移约束，解出，引入位移约束，解出 未知的结点位移分量未知的结点位移分量 d d ，并进一步求得结构的应力。，并进一步求得结构的应力。6-8 6-8 解题的具体步骤解题的具体步骤 单元的划分单元的划分应用有限单元法解题的步骤：应用有限单元法解题的步骤：(1)将计算对象划分成许多三角形单元，并按一定规律将所有的结点和单元将计算对象划分成许多三角形单元，并按一定规律将所有的结点和单元 分别编上号码；分别编上号码；(2)选定一个直角坐标系，进行单元分析，得到各单元的单元劲度矩阵选定一个直角坐标系，进行单元分析，得到各单元的单元劲度矩阵 k e 和结点荷载列阵和结点荷载列阵 FL e ；(3) 进行整体分析，组集整体劲度矩阵进行整体分析，组集整体劲度矩阵 K 和整体结点荷载列阵和整体结点荷载列阵 FL 。