《有限元分析与应用》

第1讲 引论 力学的分类

《有限元分析与应用》，曾攀，清华大学出版社

• 质点描述轨道——质点力学
• 刚体描述姿态——刚体力学
• 变形体描述姿态的耦合影响——弹性力学
• 复杂形状变形——弹塑性力学

有限元分析的力学基础是弹性力学，方程求解的数学原理是加权残值法和泛函极值原理，实现的方法是数值化离散技术，最终的载体是有限元分析软件。

注：为书写和理解统一，向量在本文中有时又称为矩阵

变形体力学

三个方面：力的平衡、变形状态、材料行为
三大变量：应力\(\sigma (x)\)、应变\(\epsilon\)、位移
三大方程：平衡方程，物理方程，几何方程
边界：位移边界，力边界

应力-应变曲线：哑铃状样品按标准执行

一维线性力学模型——三大方程

• 平衡方程：\(\sigma (x)=\frac{F}{A}\)
• 物理方程：\(\sigma (x)=E\epsilon\)
• 几何方程：\(\epsilon (x)=\frac{u}{l}\)

微分方程求解

一维杆模型

\(p\)为力密度，是一个常量，拉力\(F=\int_{0}^{L} p dx\)

由\(\sigma = E\epsilon\)可知，任意长度微元\(dx\)内都有，\(\frac {\int_{x}^{L} p dx}{A}=E\frac{du}{dx}\)

故
微分方程：\(\frac{d^2u}{dx^2}+\frac{p}{A\cdot E}=0\)
边界条件（左端固定，右端自由）：\(u|_{x=0}=0，\frac{du}{dx}|_{x=L}=0\)

方程求解

解析法：\(u(x)=-\frac{1}{2}\frac{p}{AE}x^2+\frac{p}{AE}Lx\)

差分法：

采用向前差分\(\frac{u_{i-1}-2u_i+u_{i+1}}{\Delta l^2}=\frac{d^2u}{dx^2}\)

1. 划分五个网格单元，中间四个节点

\(\left\{\begin{matrix}u_0-2u_1+u_2+\frac{1}{5^2}\frac{PL^2}{A\cdot E} =0,节点1\\ u_1-2u_2+u_3+\frac{1}{5^2}\frac{PL^2}{A\cdot E} =0,节点2\\ u_2-2u_3+u_4+\frac{1}{5^2}\frac{PL^2}{A\cdot E} =0,节点3\\u_3-2u_4+u_5+\frac{1}{5^2}\frac{PL^2}{A\cdot E} =0,节点4\end{matrix}\right.\)

边界：\(u|_{x=0}=u_0=0，u_5=u_4\)，代入方程中

\(\begin{pmatrix} -2 & 1 & 0 & 0\\ 1 & -2 & 1 & 0\\ 0 & 1 & -2 & 1\\ 0 & 0 & 1 & -1\end{pmatrix}\begin{pmatrix} u_1\\u_2 \\u_3 \\u_4\end{pmatrix}+\frac{1}{5^2} \frac{PL^2}{AE}\begin{pmatrix}1 \\1 \\1 \\1\end{pmatrix} =0\)

矩阵方程解为：

\(\begin{pmatrix} u_1 & u_2 & u_3 & u_4\end{pmatrix}^T=\frac{PL^2}{AE} \begin{pmatrix} 0.25 & 0.4375 & 0.5623 & 0.625\end{pmatrix}^T\)

2. 划分n个网格单元，中间n-1个节点

矩阵方程：

\(\begin{pmatrix} -2 & 1 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & 0\\ 0 & 0 & \ddots & \ddots & \ddots & 0\\ 0 & 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 1 & 1\end{pmatrix}\begin{pmatrix} u_1\\u_2 \\u_3 \\\vdots \\u_{n-2} \\u_{n-1}\end{pmatrix}+\frac{1}{n^2}\frac{PL^2}{AE} \begin{pmatrix} 1\\1 \\1 \\\vdots \\1 \\1\end{pmatrix}=0\)

边界：\(u|_{x=0}=u_0=0，u_n=u_{n-1}\)，代入方程中

差分法和解析法结果对比

分段数n	L/5	2L/5	3L/5	4L/5	L
5	0.250	0.438	0.563	0.625	0.625
50	0.185	0.329	0.431	0.491	0.510
100	0.183	0.324	0.425	0.486	0.505
500	0.180	0.321	0.421	0.481	0.501
解析解	0.180	0.320	0.420	0.480	0.500
网格划分越多，差分解越接近解析解

试函数法：伽辽金法-加权余量法

• 假定一个含有位置系数的函数族\(\bar u(x)\)（基函数为幂函数或者三角函数），该函数满足边界条件
• 将试函数代入控制方程中，使得残差函数最小来确定试函数的系数
• 残差函数（余量）通过权函数（函数族）得到最小值

已知

微分方程：\(\frac{d^2u}{dx^2}+\frac{p}{A\cdot E}=0\)
边界条件（左端固定，右端自由）：\(u|_{x=0}=0，\frac{du}{dx}|_{x=L}=0\)

求解

\(\bar u(x)=c_1\phi _1(x)=c_1 \sin \frac{\pi x}{2L}\)，显然满足边界条件

将试函数代入控制微分方程中，得到残差函数\(R(x)\)：
\(R(x)=\frac{d^2 \bar u _1}{dx^2}+\frac{p}{AE}=-C_1(\frac{\pi}{2L})^2 \sin \frac{\pi x}{2L}+\frac{p}{AE} \ne 0\)

使用族函数加权，求和（积分）：
\(I=\int_{0}^{L} \phi _1(x)R(x)dx=\int_{0}^{L} c_1 \sin \frac{\pi x}{2L} R(x)dx=0\)

确定待定系数\(c_1=\frac{ 16}{\pi ^3}\frac{PL^2}{AE}\)

代回试函数，\(\bar u(x)=c_1\phi _1(x)=\frac{ 16}{\pi ^3}\frac{PL^2}{AE} \sin \frac{\pi x}{2L}\)

也可以使用多待定系数

试函数	L/5	2L/5	3L/5	4L/5	L
\(\bar u(x)_1\)	0.159	0.303	0.417	0.491	0.516
\(\bar u(x)_2\)	0.175	0.321	0.423	0.480	0.497
解析解	0.180	0.320	0.420	0.480	0.500

微分方程求解方法对比

	解析解	差法法	试函数法
求解过程	解微分方程	微商代替微分；线性方程组；求解	满足边界的试函数；代回控制方程；余量最小；线性方程组
技术关键	满足全场条件的解函数		全场试函数只需满足一定边界条件
难易程度	很难	较简单	简单
求解精度	极高	高，计算量大	较高，计算量小
方法规范性	不规范，技巧不统一	规范	试函数确定不规范，后续规范
解题范围	一定范围	涵盖复杂领域	范围大

函数逼近

以为函数\(f(x)\)，\(x\in [x_0,x_L]\)

基于全域展开
傅里叶级数展开，\(f(x)=\sum _{i=0}^{n} c_i \phi _i(x\in[x_0,x_L])\)，其中\(\phi _i\)为基底函数，\(c_i\)为系数

基于子域展开
采用线性函数在子域\([x_i,x_{i+1}\)上分段展开，\(f(x)=\sum _{i=0}^{n} \{a_i+b_i x (x\in[x_i,x_{i+1})\}\)，其中\(a_i,b_i\)为系数，基底函数为一次函数

	全域展开	子域展开
优点	基底函数高次连续，容易逼近高精度	基函数简单；原始微分方程可转化为线性方程组
缺点	基底函数复杂，多维较难	基函数连续性阶次低，分段区间多

复杂几何域的函数表征与逼近

非规则的几何域很难直接构造全域的试函数，进行网格划分，网格片区构造可实现，然后进行拼接成全域

在网格单元上构建试函数，并建立力学方程，以处理复杂几何问题

单元——有限元中构建复杂结构的基本组成

• 1D：梁，杆，弹簧
• 2D：四边形，三角形，曲边六边形，八节点四边形（二阶四边形）
• 3D：八节点六面体，多节点曲面体
• 板壳：结构壳单元，热学壳单元

有限元方法和软件发展

两条路线找试函数，最后发展成有限元

• 工程：1870年变分法，1940年相似结构置换，1955年直接连续单元（矩阵力学）
• 数学：1795和1915年加权残值法，有限差值法，可变的有限差值法

有限元研究开拓者

• RICHARD COURANT，美国数学家
• JOHNARGYRIS，德国人，计算机有限元计算
• OLGIERD CECIL ZIENKIEWIZ，英国人，第一本有限元领域专著《有限元方法》

主要软件

• 1966：nastran
• 1969：Marc
• 1970：Ansys
• 1975：ADINA
• 1978：DYNA
• 1979：Abaqus
• 1982：cosmos
• 1984：ALGOR

第2讲 直接刚度法—弹簧力学分析

弹簧的力学分析原理

遵循胡克定律：力正比于弹簧位移

对于弹簧网格单元

那么
\(\left\{\begin{matrix} k(u_2-u_1)=F_2\\F_1+F_2=0\end{matrix}\right.\)

转换为场方程
\(\left\{\begin{matrix}k(u_1-u_2)=F_1\\ k(u_2-u_1)=F_2\end{matrix}\right.\)

矩阵形式，\(K U=F\)——平衡方程
\(\begin{bmatrix} k&-k \\-k &k\end{bmatrix}\begin{bmatrix} u_1\\u_2\end{bmatrix}=\begin{bmatrix} F_1\\F_2\end{bmatrix}\)

刚度矩阵组装

\(\begin{bmatrix} k_1&-k_1 &0\\-k_1 &k_1+k_2&-k_2\\0&-k_2&k_2\end{bmatrix}\begin{bmatrix} u_1\\u_2\\u_3\end{bmatrix}=\begin{bmatrix} F_{R1}\\F_2\\F_{R3}\end{bmatrix}\)

边界条件：\(u_1=0,u_3=0\)

矩阵求解：\(u_2=\frac{F_2}{k_1+k_2}\)，\(F_{R1}=-k_1u_2=-\frac{k_1F_2}{k_1+k_2}\)，\(F_{R3}=-k_2u_2=-\frac{k_2F_2}{k_1+k_2}\)

弹簧单元与杆单元的比较

弹簧：\(F=k \cdot \delta\)
拉杆：\(\sigma = E\cdot \epsilon\)，即\(F/A = E\cdot \delta/l\)

故：\(k=\frac{EA}{l}\)，拉杆与弹簧是等价的，刚度系数

杆单元的坐标变换

局部坐标转化为全局坐标，以便对单元进行集成和组装

对于平面杆，全局坐标系\(\bar x o \bar y\)，局部坐标系\(xo\)

局部坐标系中的节点位移：\(q^e=[u_1,u_2]\)

全局坐标系中的节点位移：\(\bar q^e=[\bar u_1,\bar v_1,\bar u_2,\bar v_2]\)

存在关系：\(u_1=\bar u_1 \cos \alpha +\bar v_1 \sin \alpha\)，\(u_2=\bar u_2 \cos \alpha +\bar v_2 \sin \alpha\)

\(q^e=[u_1,u_2]=\begin{bmatrix}\cos \alpha & \sin \alpha & 0 & 0\\ 0 & 0 & \cos \alpha & \sin \alpha \end{bmatrix}\begin{bmatrix} \bar u_1\\ \bar v_1\\\bar u_2 \\\bar v_2\end{bmatrix}=T^e \cdot \bar q ^e\)

代入平衡方程：\(K^eq^e=F^e\)

得到：\(K^e(T^e \cdot \bar q ^e)=F^e\)，两边左乘旋转矩阵的转置\(T^{eT}\)

即：\(\bar K^e\cdot \bar q^e=\bar F^e\)，其中\(\bar K^e=T^{eT} K^e T^e\)和\(\bar F^e=T^{eT}F^e\)表示全局坐标系下的单元刚度矩阵和节点力矩阵

同理，空间杆单元的坐标变换

四杆结构的实例分析

已知四杆相同材质和横截面积\(A=100mm^2\)，弹性模量\(E=29.5 \times 10^4 N/mm^2\)

1. 结构的离散化

节点标记及坐标

Node	x	y
1	0	0
2	400	0
3	400	300
4	0	300

单元编号以及对应的节点、长度和轴线的方向余弦

Element	start	end	l	\(n_x\)	\(n_y\)
1	1	2	400	1	0
2	3	2	300	0	-1
3	1	3	500	0.8	0.6
4	4	3	400	1	0

2. 各个单元的矩阵描述

全局坐标系下各个单元的刚度矩阵

3. 刚度矩阵的组装

对应位移的矩阵元素进行累加，组成刚度矩阵和节点载荷

刚度矩阵：\(K=K^{(1)}+K^{(2)}+K^{(3)}+K^{(4)}\)
节点位移：\(q=\begin{bmatrix} u_1 & v_1 & u_2 & v_2 & u_3 & v_3 & u_4 & v_4\end{bmatrix}^T\)
节点力：\(P=R+F=\begin{bmatrix} R_{x1} & R_{y1} & 2\times 10^4 & R_{y2} & 0 & -2.5\times 10^4 & R_{x4} &R_{y4}\end{bmatrix}^T\)

其中，\((R_{x1},R_{y1})\)为节点1处沿x和y方向的支反力，\(R_{y2}\)为节点2处沿y方向的支反力，\((R_{x4},R_{y4})\)为节点4处沿x和y方向的支反力，

组装成整体的刚度方程

4. 边界条件的处理和刚度方程求解

边界条件BC：\(u_1=v_1=v_2=u_4=v_4=0\)，代入刚度方程，化简

求解：\(\begin{bmatrix} u_2\\u_3 \\v_3\end{bmatrix}=\begin{bmatrix} 0.2712\\0.0565 \\-0.2225\end{bmatrix}mm\)

故\(q=\begin{bmatrix} 0 & 0 & 0.2712 & 0 & 0.0565 & -0.225 & 0 & 0\end{bmatrix}^T mm\)

5. 支反力的计算

四杆结构的ANSYS实例分析

软件仿真流程

• 分析类型和单元类型
• 前处理，材料，几何模型，网格划分
• 边界条件，求解
• 检查结果并处理

使用杆单元link建立模型

第3讲 针对复杂几何形状变形体的力学描述(1)

力学描述的基本思路及关于变形体材料的基本假设

材料力学

• 拉伸：杆
• 扭转：扭杆
• 弯曲：梁

弹性体基本假设5条：

• 物体内的物质连续性continuity假设，可使用连续函数描述
• 物体内的物质均匀性homogeneity假设，各个位置具有相同特性
• 物体内的物质特性各向同性isotropy假设，各个方向具有相同的力学特性
• 线弹性 linear elasticity假设，变形遵循胡克定律，去除外力后，物体可恢复
• 小变形small deformation假设，物体变形量远远小于物体几何尺寸，即在建立方程时，可以忽略高阶小量（二阶以上）

指标记法

坐标系为右手

矢量的分量形式：\(F=[F_1,F_2,F_3]=F_i,(i=1,2,3)\)

• 自由指标，每项中只出现一次的下表
  • \(\sigma _{ij}\)中的\(i,j\)在(1,2,3)中只出现一次
• 哑指标，针对乘法时的指标默认处理
  • \(\sigma _{ij}x_j=b_i\)中\(j\)为哑指标，可以在(1,2,3)变化
• Einstein求和约定，哑指标意味着求和
  • \(\sum _{j=1}^3a_{ij}x_j=b_j,(i=1,2,3)\)可以简化为\(a_{ij}x_j=b_j,(i=1,2,3)\);等价于三个线性方程的组合
• Voigt标记，将高阶自由指标的张量遵循移动规则改写成低阶张量形式
  • 对称张量，仅写一半；沿着主对角后逆时针排序元素
  • 将3*3矩阵按照角标11-22-33-23-13-12将一半元素排列成向量

三大变量及三大方程

单元的三大变量：位移，应力，应变
单元的三大方程：平衡方程，几何方程，物理方程
边界：力边界和位移边界

法应力\(\sigma\)和切(剪)应力\(\tau\)

平面问题的平衡方程构建

• 约定：正面正向为正，负面负向为正

• 有四个侧面，在平衡方程中，应考虑合力的平衡：

  • 沿方向所有合力的平衡
  • 沿y方向所有合力的平衡
  • 所有合力关于任一点的力矩平衡

• 应力在经过dx或dy变化后的位置上有增量表达

• 应力在各个侧面上为均匀分布

• 对于均匀厚度\(t\)的平板，微元\(dxdy\_t\)

• 面的法应力为\(\sigma\)，正方向为单元内部垂直指向面的外部

• 剪应力为\(\tau\)，方向垂直于法应力

• 体积力为\(b\)表示单位体积内的力载荷

1. 力学变量的增量计算

泰勒展开
\(\sigma _{xx}(x+dx,y)=\sigma _{xx}(x,y)+\frac{\partial \sigma _{xx}(x,y)}{\partial x} dx+\frac{\partial ^2 \sigma _{xx}(x,y)}{2\partial x^2} dx^2+…\)

略去二阶以上微量
\(\sigma _{xx}(x+dx,y)=\sigma _{xx}(x,y)+\frac{\partial \sigma _{xx}(x,y)}{\partial x} dx\)

2. 合力的平衡

\(\sum F_x =0\)，\(\sum F_y =0\)，\(\sum M_0 =0\)，

x方向的力：左右面上的法应力*面积+上下面的剪应力*面积+体积力*微元体积
\(\sigma _{xx}(x+dx,y) \cdot dy \cdot t-\sigma _{xx}(x,y)\cdot dy \cdot t+\tau _{xy}(x,y+dy) \cdot dx \cdot t-\tau _{xy}(x,y) \cdot dx \cdot t+\bar b_xdxdyt=0\)
位置增量（x+dx或者y+dy）上的应力（即\(\sigma _{xx}(x+dx,y)\)和\(\tau _{xy}(x,y+dy)\)）使用一阶泰勒展开替换

\(\frac{\partial \sigma _{xx}(x,y)}{\partial x} dxdyt+\frac{\partial \tau _{xy}(x,y)}{\partial y} dydxt+\bar b_xdxdyt=0\)

等式两边除以微元体积\(dxdyt\)，函数对xy连续且可微，雅可比行列式为

\(\frac{\partial \sigma _{xx}(x,y)}{\partial x} +\frac{\partial \tau _{xy}(x,y)}{\partial y} +\bar b_x=0\)，简记为\(\frac{\partial \sigma _{xx}}{\partial x} +\frac{\partial \tau _{xy}}{\partial y} +\bar b_x=0\)

同理，y方向上的合力平衡

\(\frac{\partial \sigma _{yy}}{\partial y} +\frac{\partial \tau _{yx}}{\partial x} +\bar b_y=0\)

最后，力矩平衡，即力关于微元的几何中心形成力矩平衡，对于几何中心，上下面剪应力及其偏导形成的力矩=左右面剪应力及其偏导形成的力矩；（体积力和法应力的力臂为0）
\(\tau _{xy}(x,y)dxt \cdot \frac{y}{2} +\tau _{xy}(x,y+dy)dxt \cdot \frac{y}{2} =\tau _{yx}(x+dx,y)dyt \cdot \frac{x}{2} +\tau _{yx}(x,y)dyt \cdot \frac{x}{2}\)
同样位置增量的剪应力使用一阶泰勒展开替换，则有：

\(\tau _{xy}=\tau _{yx}\)，即剪应力互等定理

汇总三个力平衡方程

\(\left\{\begin{matrix} \frac{\partial \sigma _{xx}}{\partial x} +\frac{\partial \tau _{xy}}{\partial y} +\bar b_x=0\\ \frac{\partial \sigma _{yy}}{\partial y} +\frac{\partial \tau _{yx}}{\partial x} +\bar b_y=0\\\tau _{xy}=\tau _{yx}\end{matrix}\right.\)

剪应力互等，故

\(\left\{\begin{matrix} \frac{\partial \sigma _{xx}}{\partial x} +\frac{\partial \tau _{xy}}{\partial y} +\bar b_x=0\\ \frac{\partial \sigma _{yy}}{\partial y} +\frac{\partial \tau _{xy}}{\partial x} +\bar b_y=0\end{matrix}\right.\)

哑指标形式：\(\sigma _{ij,j}+\bar b_i=0\)，其中\(\sigma _{ij,j}\)的下角标\((ij,j)\)表示物理量\(\sigma _{ij}\)对\(j\)方向求偏导，i=x或者y

平面问题的几何方程构建

平面变形量的定义需要考感两个方面

• 沿各个方向上的长度相对变化量（应变）
• 夹角的变化量

位移场\(u(x,y),v(x,y)\)分别表示任意位置的沿x和y方向的位移量

1. 沿各个方向上的应变（即单位长度的变化量）

忽略\(P'A'\)的小角度变化，那么

\(P'A'=dx+(u+\frac{\partial u}{\partial x}dx)-u=dx+\frac{\partial u}{\partial x}dx\)
\(PA=dx\)
x方向上的应变\(\epsilon _{xx} = \frac{P'A'-PA}{PA}=\frac{\partial u}{\partial x}\)
同理，y方向上的应变\(\epsilon _{yy} = \frac{P'B'-PB}{PB}=\frac{\partial v}{\partial y}\)

2. 夹角变化量

\(P'A'\)与\(PA\)的夹角\(\alpha = \tan \alpha =\frac {v+\frac{\partial v}{\partial x}dx-v}{dx}=\frac{\partial v}{\partial x}\)
\(P'B'\)与\(PB\)的夹角\(\beta = \tan \beta =\frac {u+\frac{\partial u}{\partial y}dy-u}{dy}=\frac{\partial u}{\partial y}\)
定义夹角的总变化为\(\gamma_{xy}= \alpha+\beta=\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\)

汇总平面问题的几何方程

\(\left\{\begin{matrix} \epsilon _{xx} =\frac{\partial u}{\partial x} \\ \epsilon _{yy} =\frac{\partial v}{\partial y}\\\gamma_{xy}= \frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\end{matrix}\right.\)

指标形式：\(\epsilon _{ij}=\frac{1}{2}(u_{i,j}+u_{j.i})\)

其中，\(\epsilon _{ij}=\frac{1}{2}\gamma_{i,j}(i\ne j)\)

就平面问题

• 如果已知2个位移分量u和v，可以通过式唯一求出3个应变分量\(\epsilon _{xx}，\epsilon _{xy}，\epsilon _{yy}\)
• 但如果是一个反问题，即已知3个应变分量是\(\epsilon _{xx}，\epsilon _{xy}，\epsilon _{yy}\)，就不一定能够惟一求出2个位移分量u和v，除非这3个应变分量满足一定的关系--变形协调条件（compatibilitycondition）\(\frac{\partial ^2 \epsilon _{xx}}{\partial y^2} +\frac{\partial ^2 \epsilon _{yy}}{\partial x^2}=\frac{\partial ^2 \gamma _{xy}}{\partial x\partial y}\)
• 变形协调条件的物理意义是，材料在变形过程中应是整体连续的，不应出现“撕裂”和“重叠”现象

平面问题的物理方程构建

定义x方向为主（拉伸）方向上的伸长：\(\epsilon _x= \frac{\sigma _x} {E}\)

泊松比定义为垂直拉伸方向的缩短与拉伸方向的伸长之比：\(\mu = -\frac{\epsilon_y}{\epsilon_x}\)

对于平面的微元\(dxdy\_t\)，法应力和剪应力可分解为拉伸作用和剪切作用

x方向的主拉伸\(\sigma _{xx}\)会产生泊松效应

\(\epsilon _{xx}= \frac{\sigma _{xx}} {E},\epsilon _{yy}= -\frac{\mu \sigma _{xx}} {E}\)

同理，y方向的主拉伸\(\sigma _{xx}\)也会产生泊松效应，\(\epsilon _{xx}= -\frac{\mu \sigma _{yy}} {E},\epsilon _{yy}= \frac{\sigma _{yy}} {E}\)

剪切力只产生角度变化，\(\gamma _{xy}= \frac{\tau _{xy}} {G}\)

组装成本构物理方程——本构模型，又称材料的力学 本构方程 ，或材料的应力-应变模型，是描述材料的力学特性 （应力-应变-强度-时间关系）的数学表达式

\(\left\{\begin{matrix} \epsilon _{xx} =\frac{1}{E}(\sigma _{xx}-\mu \sigma _{yy}) \\ \epsilon _{yy} =\frac{1}{E}(\sigma _{yy}-\mu \sigma _{xx})\\\gamma_{xy}= \frac{1}{G}\tau _{xy}\end{matrix}\right.\)

或者逆形式

\(\left\{\begin{matrix} \sigma _{xx} =\frac{E}{1-\mu ^2}(\epsilon _{xx}+\mu \epsilon _{yy}) \\ \sigma _{yy} =\frac{E}{1-\mu ^2}(\epsilon _{yy}+\mu \epsilon _{xx})\\\tau_{xy}= G\gamma _{xy}\end{matrix}\right.\)

其中\(G=\frac{E}{2(1+\mu)}\)

矩阵形式

\(\begin{bmatrix} \sigma _{xx} \\\sigma _{yy} \\\tau _{xy} \end{bmatrix}=\begin{bmatrix} \frac{E}{1-\mu ^2} & \frac{E\mu}{1-\mu ^2} & 0\\ \frac{E\mu}{1-\mu ^2} & \frac{E}{1-\mu ^2}& 0\\0 & 0 &G\end{bmatrix}\begin{bmatrix} \epsilon _{xx} \\\epsilon _{yy} \\\gamma _{xy} \end{bmatrix}\)

两类边界条件

BC（boundary condition）：外部描述，包括位移方面和力平衡方面的边界条件，即变形体的几何空间\(\Omega\)其表面\(\partial \Omega\)被位移边界\(S_u\)和力边界\(S_p\)完全不重叠地包围

• \(\Omega=S_u+S_p\)
• \(S_u\cap S_p=0\)

给定\(\bar u_i,\bar p_i\)

• \(u_i=\bar u_i\)，On \(S_u\)
• On \(S_p\)，\(n_x=dy/ds,n_y=dx/ds\)
  • \(\sigma _{xx} \cdot n_x+\tau _{xy} \cdot n_y=\bar p_x\)
  • \(\sigma _{yy} \cdot n_y+\tau _{yx} \cdot n_x=\bar p_y\)
  • \(\tau _{xy}=\tau _{yx}\)

※ 第4讲 针对复杂几何形状变形体的力学描述(2) ※

三大变量

• 位移：\(u(x,y),v(x,y)\)
• 应力：\(\sigma _{xx}(x,y),\sigma _{yy}(x,y),\tau _{xy}(x,y)\)
• 应变：\(\epsilon _{xx}(x,y),\epsilon _{yy}(x,y),\gamma _{xy}(x,y)\)

三大方程

1. 平衡方程

• \(\left\{\begin{matrix} \frac{\partial \sigma _{xx}}{\partial x} +\frac{\partial \tau _{xy}}{\partial y} +\bar b_x=0\\ \frac{\partial \sigma _{yy}}{\partial y} +\frac{\partial \tau _{xy}}{\partial x} +\bar b_y=0\end{matrix}\right.\)

2. 几何方程

• \(\left\{\begin{matrix} \epsilon _{xx} =\frac{\partial u}{\partial x} \\ \epsilon _{yy} =\frac{\partial v}{\partial y}\\\gamma_{xy}= \frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\end{matrix}\right.\)

3. 物理方程

• \(\left\{\begin{matrix} \epsilon _{xx} =\frac{1}{E}(\sigma _{xx}-\mu \sigma _{yy}) \\ \epsilon _{yy} =\frac{1}{E}(\sigma _{yy}-\mu \sigma _{xx})\\\gamma_{xy}= \frac{1}{G}\tau _{xy}\end{matrix}\right.\)

两类边界条件

• On \(S_u \left\{\begin{matrix} u=\bar u\\v=\bar v\end{matrix}\right.\)
• On \(S_p \left\{\begin{matrix}\sigma _{xx} \cdot n_x+\tau _{xy} \cdot n_y=\bar p_x\\\sigma _{yy} \cdot n_y+\tau _{xy} \cdot n_x=\bar p_y\end{matrix}\right.\)

几种特殊情况的讨论

平面应力：近似认为\(\sigma _{zz}=0\)

平面应变：近似认为\(\epsilon _{zz}=0\)

平面应力（平衡方程）和平面应变（几何方程）的变量和方程完全相同，二者通过物理方程转换

刚体位移：物体内不产生任何应变，对几何方程给定零位移

故平面刚体位移的表达（\(u_0,v_0,\omega _0\)均为常数）：

\(\left\{\begin{matrix} u(x,y)=-\omega _0 y+u_0\\v(x,y)=\omega _0 x+v_0\end{matrix}\right.\)

• \(u_0\)为整个物体在x方向的刚体位移量
• \(v_0\)为整个物体在y方向的刚体位移量
• \(\omega _0\)为整个物体在空间的刚体旋转角度

简单拉杆问题的完整弹性力学求解

1D弹性杆的方程解

\(\left\{\begin{matrix} \sigma _{x}(x) =\frac{P}{A}\\ \epsilon _{x}(x) =\frac{P}{EA}\\ u(x)= \frac{P}{EA}x\end{matrix}\right.\)

• 解析方法的求解过程严谨，可得到物体内各点力学变量的表达，是场变量
• 经验方法的求解过程较简单，但需要事先进行假定，往往只得到一些特定位置的力学变量表达，而只能应用于一些简单情形

平面纯弯梁的描述及求解

假设：

• 直法线假设
• 小变形

位移：\(v(x,\tilde{y} =0)\) （中性层绕度）
应力： \(\sigma_x(x,y)\) （其它应力分量很小，不考虑），该变量对应于梁截面上的弯矩\(M\)
应变：\(\epsilon_x(x,y)\)（满足直线假定）

物理方程：\(\sigma _x =E \epsilon _x\)

※空间弹性问题的完整描述※

对于空间中的微元体\(dxdydz\)，维度扩展

基本变量

	1D	2D	3D
位移分量	\(u(x)\)	\(u(x,y),v(x,y)\)	\(u(x,y,z),v(x,y,z),w(x,y,z)\)
应变分量	\(\epsilon _{xx}(x)\)	\(\epsilon _{xx}(x,y),\epsilon _{yy}(x,y),\gamma _{xy}(x,y)\)	\(\epsilon _{xx}(x,y,z),\epsilon _{yy}(x,y,z),\epsilon _{xy}(x,y,z)\) \(\gamma _{xy}(x,y,z),\gamma _{yz}(x,y,z),\gamma _{xz}(x,y,z)\)
应力分量	\(\sigma _{xx}(x)\)	\(\sigma _{xx}(x,y),\sigma _{yy}(x,y),\tau _{xy}(x,y)\)	\(\sigma _{xx}(x,y,z),\sigma _{yy}(x,y,z),\sigma _{xy}(x,y,z)\) \(\tau _{xy}(x,y,z),\tau _{yz}(x,y,z),\tau _{xz}(x,y,z)\)
平衡方程	\(\frac{\partial \sigma _{xx}}{\partial x} =0\)	\(\left\{\begin{matrix} \frac{\partial \sigma _{xx}}{\partial x} +\frac{\partial \tau _{xy}}{\partial y} +\bar b_x=0\\ \frac{\partial \sigma _{yy}}{\partial y} +\frac{\partial \tau _{xy}}{\partial x} +\bar b_y=0\end{matrix}\right.\)	\(\left\{\begin{matrix} \frac{\partial \sigma _{xx}}{\partial x} +\frac{\partial \tau _{xy}}{\partial y}+\frac{\partial \tau _{xz}}{\partial z} +\bar b_x=0\\ \frac{\partial \tau _{xy}}{\partial x}+\frac{\partial \sigma _{yy}}{\partial y} +\frac{\partial \tau _{yz}}{\partial z} +\bar b_y=0\\ \frac{\partial \tau _{xz}}{\partial x} +\frac{\partial \tau _{yz}}{\partial y} +\frac{\partial \sigma _{zz}}{\partial z} +\bar b_z =0 \end{matrix}\right.\)
几何方程	\(\epsilon _{xx} =\frac{d u}{d x}\)	\(\left\{\begin{matrix} \epsilon _{xx} =\frac{\partial u}{\partial x} \\ \epsilon _{yy} =\frac{\partial v}{\partial y}\\\gamma_{xy}= \frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\end{matrix}\right.\)	\(\left\{\begin{matrix} \epsilon _{xx} =\frac{\partial u}{\partial x} , \epsilon _{yy} =\frac{\partial v}{\partial y}, \epsilon _{zz} =\frac{\partial w}{\partial z}\\ \gamma_{xy}= \frac{\partial v}{\partial x}+\frac{\partial u}{\partial y},\gamma_{yz}= \frac{\partial w}{\partial y}+\frac{\partial v}{\partial z},\gamma_{xz}= \frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\end{matrix}\right.\)
物理方程	\(\epsilon _{xx} =\frac{\sigma _{xx}}{E}\)	\(\left\{\begin{matrix} \epsilon _{xx} =\frac{1}{E}(\sigma _{xx}-\mu \sigma _{yy}) \\ \epsilon _{yy} =\frac{1}{E}(\sigma _{yy}-\mu \sigma _{xx})\\\gamma_{xy}= \frac{1}{G}\tau _{xy}\end{matrix}\right.\)	\(\left\{\begin{matrix} \epsilon _{xx} =\frac{1}{E}(\sigma _{xx}-\mu (\sigma _{yy}+\sigma _{zz})) \\ \epsilon _{yy} =\frac{1}{E}(\sigma _{yy}-\mu (\sigma _{xx}+\sigma _{zz}))\\ \epsilon _{zz} =\frac{1}{E}(\sigma _{zz}-\mu (\sigma _{xx}+\sigma _{yy})) \\ \gamma_{xy}= \frac{1}{G}\tau _{xy},\gamma_{yz}= \frac{1}{G}\tau _{yz},\gamma_{xz}= \frac{1}{G}\tau _{xz} \end{matrix}\right.\)
几何边界\(BC(u)\)	\(u(x)|_{x=x_0}=\bar u\)	\(u(x,y)|_{x=x_0,y=y_0}=\bar u\) \(v(x,y)|_{x=x_0,y=y_0}=\bar v\)	\(u(x,y,z)|_{x=x_0,y=y_0,z=z_0}=\bar u\) \(v(x,y,z)|_{x=x_0,y=y_0,z=z_0}=\bar v\) \(w(x,y,z)|_{x=x_0,y=y_0,z=z_0}=\bar w\)
外力边界\(BC(p)\)	\(\sigma(x)|_{x=x_0}=\bar p_x\)	\(\sigma _{xx}(x_0,y_0) n_x+\tau _{xy}x_0,y_0)n_y=\bar p_x\) \(\sigma _{yy}(x_0,y_0) n_y+\tau _{yx} (x_0,y_0) n_x=\bar p_y\)	\(\sigma _{xx}(x_0,y_0,z_0) n_x+\tau _{xy}(x_0,y_0,z_0)n_y+\tau _{xz}(x_0,y_0,z_0)n_z=\bar p_x\) \(\tau _{xy}(x_0,y_0,z_0) n_x+\sigma _{yy}(x_0,y_0,z_0)n_y+\sigma _{xz}(x_0,y_0,z_0)n_z=\bar p_y\) \(\tau _{xz}(x_0,y_0,z_0) n_x+\tau _{yz}(x_0,y_0,z_0)n_y+\sigma _{zz}(x_0,y_0,z_0)n_z=\bar p_z\)

三组剪应变相等：\(\gamma _{xy}=\gamma _{yx},\gamma _{yz}=\gamma _{zy},\gamma _{zx}=\gamma _{xz}\)
三组剪应力相等：\(\tau _{xy}=\tau _{yx},\tau _{yz}=\tau _{zy},\tau _{zx}=\tau _{xz}\)
独立变量变为：3（方向位移）+6（法向应变+剪切应变）+6（法向应力+剪切应力）

• \(x_0,y_0,z_0\)为边界上的几何坐标
• \(n_x,n_y,n_z\)为边界外法线的方向余弦
• \(\bar u,\bar v,\bar w\)为给定的对应方向上的位移值
• \(\bar p_x,\bar p_y,\bar p_z\)为给定的对应方向上的边界分布力

逆形式

关于张量的描述及理解

• 0阶张量：无自由指标的量，与坐标系选取无关，如温度、质量、能量等标量
• 阶张量：有1个自由指标的量，如坐标x，位移u等矢量
• 2阶张量：有2个自由指标的量，如应力、应变
• n阶张量：有n个自由指标的量，如四阶弹性系数张量

同一矢量在不同坐标系下的分量形式不同，但模长不变——张量的不变量

应力转换

应力莫尔圆中应力最大特征量——第一特征量，最小特征量——第二特征量

• 第一和第二特征量可以作为材料强度准则

莫尔圆

应力应变莫尔圆（Mohr's Circle）是一种图示工具，用于分析二维应力或应变状态。它可以帮助工程师和材料科学家可视化不同方向上材料内部的应力和应变情况。

莫尔圆的定义：

• 莫尔圆是一个图形，用于表示在任意方向上的应力和应变状态。通过在圆内描绘不同方向下的法向应力和切向应力，可以直观地展示材料内部的应力状态。

1. 确定基本应力：

   • 标识出直应力\((\sigma_x、\sigma_y)\)和剪应力\((\tau_{xy})\)的值。

2. 画坐标系：

   • 选择一个二维坐标系，横轴表示法向应力(\(\sigma\))，纵轴表示剪应力\((\tau)\)

3. 标记点：

   • 在坐标系中标记出应力状态下的两个点：
     • 点A：\((\sigma_x, \tau_{xy})\)
     • 点B：\((\sigma_y, -\tau_{xy})\)

4. 连接点并绘制圆：

   • 连接这两个点，然后找到圆心\((\frac{\sigma_x+\sigma_y}{2},0)\)，并根据两个点的距离绘制完整的圆
   • \(\sigma ^2-(\sigma_x+\sigma_y)\sigma+\sigma_x\sigma_y+\tau ^2- \tau _{xy}^2=0\)

5. 计算主应力：

   • 圆的交点表示材料的主应力状态，主应力（\(\sigma_1\)和 \(\sigma_2\)）位于圆的极点上
   • 主应力（\(\sigma_1\)和 \(\sigma_2\)）在横轴表示法向应力(\(\sigma\))，此时剪应力为0

理解莫尔圆

• 几何意义：莫尔圆的几何特征直观地表明了法向应力和切向应力之间的关系。圆上的每一点对应一个特定的方向上的应力状态

• 主应力：通过求出圆的交点，可以获得主应力，这些应力是材料内部的“最大”与“最小”应力状态，确保在设计和分析中考虑材料的弱点

• 切向应力变化：圆的大小和位置说明了在不同方向上剪切和法向应力的变化。有助于理解材料在不同加载条件下的表现，包括屈服和破坏

第5讲 变形体力学方程求解的试函数方法的原理

变形体力学方程求解的主要方法分类及试函数方法

1、直接针对原始方程进行求解，方法有

• 解析法(analytical method)
• 半解析法(semi-inverse method）
• 差分法(finite difference method）

2、间接针对原始方程进行求解（误差处理），方法有

• 加权残值法(weighted residual method)
  • Galerkin加权残值法
  • 残值最小二乘法（least sguares method）
• 虚功原理(principle of virtual work）
• 最小势能原理(principle of minimum potential energy)
• 变分方法(variational method)

Galerkin加权残值法

残值最小二乘法

平面弯曲梁求解的试函数方法-残值处理法

控制方程：\(EI\frac{d^4v}{dx^4}=\bar p_0\)
位移边界条件：\(v(x)|_{x=0}=0\)，\(v(x)|_{x=l}=0\)
载荷边界条件：\(v''(x)|_{x=0}=0\)，\(v''(x)|_{x=l}=0\)

1. 单系数试函数

\(\widehat{v}_1(x)=c_1 \phi(x)=c_1 \cdot \sin \frac{\pi x}{l}\)

代入控制方程得到余量函数，\(\Re =EI\frac{d^4(\widehat{v}_1(x))}{dx^4}-\bar p_0 \ne 0\)

选择基底函数作为权函数\(w= \sin \frac{\pi x}{l}\)

残差方程，\(\int_{0}^{l} w \cdot \Re dx=\int_{0}^{l} (\sin \frac{\pi x}{l})[EI\frac{d^4(\widehat{v}_1(x))}{dx^4}-\bar p_0 ]dx =0\)

求解得到\(c_1=\frac{4l^4}{\pi ^5 EI}\bar p_0\)，故试函数\(\widehat{v}_1(x)=\frac{4l^4}{\pi ^5 EI}\bar p_0 \cdot \sin \frac{\pi x}{l}\)

2. 多系数试函数

\(\widehat{v}_2(x)=c_1 \phi _1(x)+c_2 \phi _2(x)=c_1 \cdot \sin \frac{\pi x}{l}+c_2\sin \frac{3 \pi x}{l}\)

代入控制方程得到余量函数，\(\Re =EI\frac{d^4(\widehat{v}_2(x))}{dx^4}-\bar p_0 \ne 0\)

选择基底函数作为权函数\(w_1=\phi_1 =\sin \frac{\pi x}{l}\)和\(w_2=\phi _2= \sin \frac{3 \pi x}{l}\)

残差方程

\(\int_{0}^{l} w_1 \cdot \Re dx=\int_{0}^{l} (\sin \frac{\pi x}{l})[EI\frac{d^4(\widehat{v}_2(x))}{dx^4}-\bar p_0 ]dx =0\)
\(\int_{0}^{l} w_2 \cdot \Re dx=\int_{0}^{l} (\sin \frac{3 \pi x}{l})[EI\frac{d^4(\widehat{v}_2(x))}{dx^4}-\bar p_0 ]dx =0\)

求解得到\(c_1=\frac{4l^4}{\pi ^5 EI}\bar p_0,c_2=\frac{4l^4}{243\pi ^5 EI}\bar p_0\)，故试函数\(\widehat{v}_2(x)=\frac{4l^4}{\pi ^5 EI}\bar p_0 \cdot \sin \frac{\pi x}{l}+\frac{4l^4}{243 \pi ^5 EI}\bar p_0 \cdot \sin \frac{3\pi x}{l}\)

3. 残值最小二乘法

选取试函数\(\widehat{v}_2(x)=c_1 \phi _1(x)+c_2 \phi _2(x)=c_1 \cdot \sin \frac{\pi x}{l}+c_2\cdot \sin \frac{3 \pi x}{l}\)

代入控制方程得到余量函数，\(\Re =EI\frac{d^4(\widehat{v}_2(x))}{dx^4}-\bar p_0 \ne 0\)

取权函数\(w_t=1\)，定义残差平方的积分为一个误差指标

\(E_{rr}=\int_{\Omega} \Re ^2(c_1,c_2,……,c_n,\phi _1,\phi _2,……,\phi_n) d \Omega\)

此时需要确定待定系数\(c_i\)，使得定义的误差指标达到最小值

由最小而二乘法，有

\(\left\{\begin{matrix}\frac {\partial E_{rr} (c_1,c_2)} {\partial c_1} =0 \\\frac {\partial E_{rr} (c_1,c_2)} {\partial c_2} =0 \end{matrix}\right.\)

关于\(c_i\)的线性方程组，解得\(c_1=\frac{4l^4}{\pi ^5 EI}\bar p_0,c_2=\frac{4l^4}{243\pi ^5 EI}\bar p_0\)，与Galerkin法的结果相同，其原因是基底函数选取是一样的

如何降低对试函数的高阶导数要求

控制方程中存在高阶导数时，使用试函数法必须保证其高阶导数也存在，虽然正弦函数可以很好地满足，其他函数选区的话比较困难

平面弹性方程，基于位移求解的控制方程及其边界条件

控制方程

边界条件

假设试函数

其中，\(c_{ui}\)和\(c_{vi}\)是待定系数，\(\phi_{ui}(x,y)\)和\(\phi _{vi}(x,y)\)为满足在所有边界条件的基底函数

将试函数代入控制方程，使其残值在加权积分下为0，即得到Garlerkin方程

得到关于\(c_{ui}\)和\(c_{vi}\)的线性方程组，问题得解

由此可知，对试函数有

• 满足所有边界条件
• 试函数二阶导数必须存在（弯曲梁，导数4阶；一般问题，为2阶）
• 全场几何域的积分计算
• 流程十分规范
• 如何降低试函数对高阶导数的要求？

能量原理：虚功原理，最小势能原理

• 只满足位移边界条件
• 降低对试函数的导数阶次要求
• 数学原理：分部积分和变分方法

平面弯曲梁求解的虚功原理

虚位移：任意扰动的位移应满足的条件（如几何方程），称为许可位移条件，把满足许可位移条件的、任意的微小位移称为虚位移(假想的）
(virtual displacement)

虚功：某点的虚位移与其负载力的乘积为虚功(virtual work)

• 平衡状态下的体系，当作用有满足许可位移条件的虚位移时，系统上的所有的虚功总和恒为零
• 系统的力分为外力和内力（应力），引起外力虚功\(\delta W\)和内力虚功\(\delta U\)，内力虚功又称为虚位移能
• 虚应变能(virtual strain energy)，由于弹性体在变形过程中，内力是克服变形所产生的，其方向总是与变形的方向相反，所以内力虚功取负，故\(\delta W=\delta U\)；这里的虚位移进满足位移边界条件的许可位移

平面弯曲梁求解的能量原理（降低试函数的导数阶次要求）

控制方程：\(EI\frac{d^4v}{dx^4}=\bar p_0\)
位移边界条件：\(v(x)|_{x=0}=0\)，\(v(x)|_{x=l}=0\)
载荷边界条件：\(v''(x)|_{x=0}=0\)，\(v''(x)|_{x=l}=0\)

1. 虚功原理求解

选取试函数\(\widehat{v}(x)=c_1 \cdot \sin \frac{\pi x}{l}\) ，显然满足位移边界条件

对待定系数\(c_1\)进行微笑变化，则虚位移场：\(\delta \widehat{v}(x)=\delta c_1 \cdot \sin \frac{\pi x}{l}\)

虚应变能：\(\delta U= \int _{\Omega} \sigma _x \delta \epsilon _x d \Omega=\int _{0}^1 \int _{A} E \epsilon_x \cdot \delta \epsilon _x dAdx\)

梁弯曲的几何方程有：\(\epsilon _x(x) =- \widehat y \frac{d^2 v(x)}{dx^2}\)；故，
\(\delta U=\int _{0}^1 E( \int _{A} \widehat y ^2 dA )\cdot (\frac{d^2 v(x)}{dx^2}) \cdot (\frac{d^2 \delta v(x)}{dx^2})dx=\frac{EIl}{2}(\frac{\pi }{l})^4 \cdot c_1 \cdot \delta c_1\)
其中截面惯性矩\(I=\int _{A} \widehat y ^2 dA\)

简支梁的外力虚功为：\(\delta W= \int _0^1 \bar p_0 \delta \widehat vdx= \bar p_0 \cdot \widehat c_1 \cdot \int _0^1\delta \sin \frac{\pi x}{l} dx=\frac{2l\bar p_0}{\pi}\cdot \delta c_1\)

虚功原理，有\(\delta W=\delta U\)：\(\frac{2l\bar p_0}{\pi}\cdot \delta c_1=\frac{EIl}{2}(\frac{\pi }{l})^4 \cdot c_1 \cdot \delta c_1\)
化简，\((\frac{2l\bar p_0}{\pi}- \frac{EIl}{2}(\frac{\pi }{l})^4 \cdot c_1 )\cdot \delta c_1=0\)

因为\(\delta c_1=0\)具有任意性，\(c_1=\frac{4l^4}{\pi ^5 EI}\bar p_0\)，故试函数\(\widehat{v}(x)=\frac{4l^4}{\pi ^5 EI}\bar p_0 \cdot \sin \frac{\pi x}{l}\)

• 该结果与前面Galerkin方法得到的结果相同，这是因为两种方法取了相同的试函数。
• 从以上的求解过程可以看出，虚功原理与加残值法类似，只能在试函数的范围内寻找最好的解，但该结果不一定是精确的
• 这取决于事先所假定的位移模式（试函数），如果事先所假定的位移模式包含有精确解的情况，那么由虚功原理求解出的解一定是精确的

2. 最小势能原理求解

选取试函数\(\widehat{v}(x)=c_1 \cdot \sin \frac{\pi x}{l}+c_2 \cdot \sin \frac{3 \pi x}{l}\) ，显然满足位移边界条件

定义势能=应变能-外力功：\(min _{\widehat v(x) \in BC(u)}[\prod (\widehat v(x))=U-W]\)
即：\(U= \frac{1}{2} \int _{\Omega} \sigma _x \cdot \epsilon _x\cdot d \Omega\)，\(W=\int _0^l \bar p_0 \cdot \widehat v(x) dx\)

真实的位移场总是使得物体的总势能取最小值

应变能\(U= \frac{1}{2} \int _{\Omega} \sigma _x \cdot \epsilon _x\cdot d \Omega=\frac{1}{2} \int _{0}^l \sigma _x \cdot (\frac{d^2 \widehat v}{dx^2})^2 \cdot dx=\frac{EI}{2}[c_1^2(\frac{\pi}{l})^4\frac{l}{2}+c_2^2(\frac{3\pi}{l})^4\frac{l}{2}]\)

外力功\(W=\int _0^l \bar p_0 \cdot \widehat v(x) dx=\int _0^l \bar p_0 \cdot (c_1 \cdot \sin \frac{\pi x}{l}+c_2 \cdot \sin \frac{3 \pi x}{l}) dx=\bar p_0[c_1\frac{2l}{\pi}+c_2\frac{2l}{3\pi}]\)
总势能\(\prod (\widehat v(x))=U-W\)，为取得极小值

\(\frac{\partial \prod}{\partial c_i}=0\)，即\(\frac{\partial \prod}{\partial c_1}=\frac{\partial \prod}{\partial c_2}=0\)

解得\(c_1=\frac{4l^4}{\pi ^5 EI}\bar p_0,c_2=\frac{4l^4}{243\pi ^5 EI}\bar p_0\)，与Galerkin法、最小二乘法以及虚功原理的结果相同，其原因是基底函数选取是一样的

最小势能原理方法也叫Rayleigh-Ritz方法（瑞利-里茨法）

平面弯曲梁求解的最小势能原理的变分基础

证明，势能最小是否等价于控制方程和载荷边界条件

可使得势能泛函取最小值的位移边界条件的一个函数

一般弹性问题的能量原理

• 虚功原理：内外力虚功相等
• 最小势能原理

1. 虚功原理

设选取了满足位移边界条件\(BC(u)\)许可位移场\(u_i\)
则虚位移\(\delta u_i\)
有几何方程知，虚应变\(\delta \epsilon _{ij}=\frac{1}{2} (\delta u_{i,j}+\delta u_{j,i})\)
虚应变能为，\(U= \frac{1}{2} \int _{\Omega} \sigma _{ij} \cdot \epsilon _{ij}\cdot d \Omega\)，\(\delta U= \int _{\Omega} \sigma _{i,j} \delta \epsilon _{i,j} d \Omega\)
外力虚功为，\(W=\int _{\Omega} \bar b_i u _i d \Omega+\int _{S_p} \bar p_i u _i dA\)，\(\delta W=\int _{\Omega} \bar b_i \delta u _i d \Omega+\int _{S_p} \bar p_i \delta u _i dA\)，其中\(\bar b_i\)和\(\bar p_i\)为体积力和面分布力
应用虚功原理\(\delta W=\delta U\)，\(\int _{\Omega} \sigma _x \delta \epsilon _x d \Omega=\int _{\Omega} \bar b_i \delta u _i d \Omega+\int _{S_p} \bar p_i \delta u _i dA\)

2. 最小势能原理

势能\(\prod (\widehat u(x))=U-W=\frac{1}{2} \int _{\Omega} \sigma _{ij} \epsilon _{ij} d \Omega-[\int _{\Omega} \bar b_i u _i d \Omega+\int _{S_p} \bar p_i u _i dA]\)

应变能中的应变和应力各有6个变量，角标相同组合成6个应变能，有两个部分

• 正应力-正应变的变形能微元\(\Delta U_{(\sigma , \epsilon)}=\frac{1}{2} \sigma _{ii} \cdot \epsilon _{ii}\cdot d \Omega\)
• 剪应力-剪应变的变形能\(\Delta U_{(\tau , \gamma)}=\frac{1}{2} \tau _{xy} \cdot \gamma _{xy}\cdot d \Omega\)
• 出现1/2是因为假设应变微元从0线性增加到\(d \epsilon\)

两种能量法是完全等价的

第6讲 基于试函数方法的经典实现及有限元实现

基于试函数的经典方法与有限元方法

经典方法：选取全域试函数
有限元方法：子域线性试函数
二者最大的不同就是是否需要离散化

有限元思路

• 几何离散化
• 单元研究，单元上选取试函数并建立力学方程
• 集成还原，刚度矩阵组装

有限元方法中的自然离散与逼近离散

自然离散

• 基于结构中的自然连接关系划分节点，并形成离散单元
• 计算精度高，甚至可以获得精确解
• 例如杆、铰链的受力问题

逼近离散

• 对于连续体结构，需要人工划分节点，以分片（单元）连续的形式来进行逼近
• 例如骨头受力
• 影响计算精度与误差的因素
  1.节点的位置和数量
  2.计算规模与计算量
  3.选择单元的类型
  4.对几何模型的逼近程度

有限元方法中的基本步骤

• 几何域的离散化
• 单元研究
• 单元的集成（组装）
• 位移边界条件的处理（已知和未知）
  
• 支反力的计算
• 单元其它物理量的计算：应力和应变

单元研究

节点描述\(\Omega = \sum \Omega ^e\)

• 参数化几何坐标（标准化，规范化）
• 节点上的位移分量
• 节点上的力分量

单元上的描述：节点描述和三大类变量的场描述

• 位移场函数
• 应变场函数
• 应力场函数

能量原理

• 最小势能原理：单元上的应变能，外力功
• 虚功原理：单元上的虚应变能，虚外力功

获得单元方程

• 对原平衡方程的加权逼近
• 对力边界条件的加权逼近

试函数经典方法及有限元方法的比较

有限元计算

• 软件平台：标准化处理，前后处理的可视化
• 硬件平台：大规模计算，非线性方程的求解

第7讲 杆、梁结构的有限元分析

局部坐标系中的杆单元构建及MATLAB编程

1. 节点基本变量及描述

上角标\(e\)表示element

杆单元具有2个节点Node，其坐标\(x_1=0，x_2=l^e\)
节点位移（向量）矩阵，\(\mathbf{q}^e=\begin{bmatrix} u_1 & u_2\end{bmatrix}^T\)
节点力（向量）矩阵，\(\mathbf{p}^e=\begin{bmatrix} P_1 & P_2\end{bmatrix}^T\)

2. 单元上的场描述

位移场函数取多项式：\(u(x)=a_0+a_1x+a_2x^2+……\)

• 唯一性原则
• 从低阶到高阶原则

杆单元含有2个节点，位移模式取\(u(x)=a_0+a_1x\)
单元节点的BC：\(u(x)|_{x=0}=u_1,u(x)|_{x=l^e}=u_2\)
解得，\(a_0=u_1,a_1=\frac{u_2-u_1}{l^e}\)
故，位移场函数，\(u(x)=u_1+\frac{u_2-u_1}{l^e}x=(1-\frac{x}{l^e})u_1+\frac{x}{l^e}u_2=\mathbf N(x) \cdot \mathbf q^e\)
其中\(\mathbf N(x) =\begin{bmatrix} 1-\frac{x}{l^e} & \frac{x}{l^e}\end{bmatrix}\)为形状函数矩阵
节点位移\(\mathbf q^e=\begin{bmatrix} u_1 & u_2\end{bmatrix}^T\)

应变场函数\(\epsilon=\frac{du(x)}{dx}=-\frac{1}{l^e}u_1+\frac{1}{l^e}u_2=\mathbf B(x) \cdot \mathbf q^e\)
其中\(\mathbf B(x) =\frac{d}{dx}\mathbf N(x)=\begin{bmatrix} -\frac{1}{l^e} & \frac{1}{l^e}\end{bmatrix}\)为几何矩阵

应力场函数\(\sigma = E^e \epsilon (x)=E^e \cdot \mathbf B(x) \cdot \mathbf q^e=\mathbf S(x) \cdot \mathbf q^e\)
其中\(\mathbf S(x) =E^e \cdot \mathbf B(x)=\begin{bmatrix} -\frac{E^e}{l^e} & \frac{E^e}{l^e}\end{bmatrix}\)为应力矩阵

3. 最小势能原理

单元的虚应变能为，\(U= \frac{1}{2} \int _{\Omega} \sigma _{ij} \cdot \epsilon _{ij}\cdot d \Omega =\frac{1}{2} \int ^{l^e}_0 [\mathbf S(x) \cdot \mathbf q^e ]^T \cdot (\mathbf B(x) \cdot \mathbf q^e) \cdot ( A^e dx)\)
\(=\frac{1}{2} E^e A^e\cdot [\mathbf B(x) \cdot \mathbf q^e ]^T \cdot (\mathbf B(x) \cdot \mathbf q^e) \cdot ( \int ^{l^e}_0 dx)\)
\(=\frac{1}{2} E^e A^e l^e\cdot \mathbf q^{eT}\cdot \mathbf B(x) ^T \cdot \mathbf B(x) \cdot \mathbf q^e\)
\(=\frac{1}{2} \frac{E^e A^e} {l^e}\cdot \mathbf q^{eT}\cdot\begin{bmatrix} 1 & -1\\-1 &1\end{bmatrix}\cdot \mathbf q^e\)

记单元刚度矩阵为\(\mathbf K^e=\frac{E^e A^e} {l^e}\begin{bmatrix} 1 & -1\\-1 &1\end{bmatrix}\)
故\(U=\frac{1}{2} \mathbf q^{eT}\cdot \mathbf K^e\cdot \mathbf q^e\)

外力虚功为，\(W=\int _{L_p} \bar P_i du=P_1u_1+P_2u_2=\mathbf P^{eT}\cdot \mathbf q^e\)，其中和\(\bar P_i\)为外力

总势能\(\prod (\mathbf q^e)=U-W=\frac{1}{2} \mathbf q^{eT}\cdot \mathbf K^e\cdot \mathbf q^e-\mathbf P^{eT}\cdot \mathbf q^e\)
泛函取极值，对\(\mathbf q^e\)取导数，\(\frac{\partial \prod (\mathbf q^e)}{\partial \mathbf q^e}=0\)
根据矩阵导数\(\nabla(Ax-b)=A^T\)和\(\nabla(Ax-b)^T(Ax-b)=2A^T(Ax-b)\)可知，\(\nabla(\frac{1}{2} \mathbf q^{eT}\cdot \mathbf K^e\cdot \mathbf q^e-\mathbf P^{eT}\cdot \mathbf q^e)=\frac{1}{2}\cdot2E^T \cdot \mathbf K^e \cdot(E\mathbf q^e)-(\mathbf P^{eT})^T=\mathbf K^e \cdot \mathbf q^e-\mathbf P^{e}=0\)
这就是杆单元的刚度方程

4. 虚功原理

节点位移\(\mathbf q^e=\begin{bmatrix} u_1 & u_2\end{bmatrix}^T\)
节点虚位移\(\delta \mathbf q^e=\begin{bmatrix} \delta u_1 & \delta u_2\end{bmatrix}^T\)
单元的虚应变能为，\(\delta U=\int _{\Omega} \sigma _x \delta \epsilon _x d \Omega= \int _{\Omega}\mathbf S(x) \cdot \mathbf q^e \cdot \mathbf B(x) \cdot \delta \mathbf q^e d \Omega\)
\(=\int ^{l^e}_0 [\mathbf S(x) \cdot \mathbf q^e ]^T \cdot (\mathbf B(x) \cdot \delta \mathbf q^e) \cdot ( A^e dx)\)
\(= \mathbf q^{eT}\cdot \mathbf K^e\cdot \delta \mathbf q^e\)

外力虚功为，\(\delta W=\bar P_i \delta u_i=P_1\cdot \delta u_1+P_2\cdot \delta u_2=\mathbf P^{eT}\cdot \delta \mathbf q^e\)，其中和\(\bar P_i\)为外力

虚功原理，\(\delta U=\delta W\)
所以，\(\mathbf q^{eT}\cdot \mathbf K^e\cdot \delta \mathbf q^e=\mathbf P^{eT}\cdot \delta \mathbf q^e\)
即，\((\mathbf q^{eT}\cdot \mathbf K^e-\mathbf P^{eT})\cdot \delta \mathbf q^e=0\)
虚位移\(\delta \mathbf q^e\)具有任意性，故\(\mathbf q^{eT}\cdot \mathbf K^e-\mathbf P^{eT}=0\)，转置后就是单元的刚度矩阵，与最小势能原理是一致的

5. 1D杆转换为平面或者空间杆，即杆单元的坐标变换

局部坐标系中的节点位移\(\mathbf q^e=\begin{bmatrix} u_1 & u_2\end{bmatrix}^T\)
2D全局坐标系中的节点位移\(\bar {\mathbf q}^e=\begin{bmatrix} \bar u_1 & \bar v_1 &\bar u_2 &\bar v_2 \end{bmatrix}^T\)
其中，\(u_1=\bar u_1 \cos \alpha+\bar v_1 \sin \alpha\)，\(u_2=\bar u_2 \cos \alpha+\bar v_2 \sin \alpha\)
转换成矩阵形式，\(\mathbf q^e=\mathbf T^e \cdot \bar {\mathbf q}^e\)
其中2D的旋转矩阵\(\mathbf T^e =\begin{bmatrix} \cos \alpha & \sin \alpha & 0 & 0\\ 0 & 0 & \cos \alpha &\sin \alpha\end{bmatrix}\)
整体坐标系下的总势能\(\prod (\mathbf q^e)=U^e-W^e=\frac{1}{2} \mathbf q^{eT}\cdot \mathbf K^e\cdot \mathbf q^e-\mathbf P^{eT}\cdot \mathbf q^e= \frac{1}{2} \bar {\mathbf q}^{eT}(\mathbf T^{eT}\mathbf K^{e} \mathbf T^{e}) {\mathbf q}^{e}-({\mathbf T}^{eT}{\mathbf P}^{e})^T\bar {\mathbf q}^{e}\)
\(=\frac{1}{2} \bar {\mathbf q}^{eT} \bar{ \mathbf K^{e} } \bar {\mathbf q}^{e}-\bar {\mathbf P}^{eT}\bar {\mathbf q}^{e}\)

其中整体刚度矩阵\(\bar {\mathbf K}^{eT}=\mathbf T^{eT}\mathbf K^{e} \mathbf T^{e}\)，整体节点力矩阵\(\bar {\mathbf P}^{e}={\mathbf T}^{eT}{\mathbf P}^{e}\)

同理，3D的旋转矩阵\(\mathbf T^e=\begin{bmatrix} \cos(x,\bar x) & \cos(x,\bar y) & \cos(x,\bar z) & 0& 0 & 0\\ 0& 0 & 0 & \cos(x,\bar x) & \cos(x,\bar y) &\cos(x,\bar z) \end{bmatrix}\)
\(\cos (x,\bar x)\)表示局部坐标系\(x\)轴与全局空间坐标系\(\bar x\)之间的方向余弦

1D杆单元分析的matlab程序

设计四个函数

• 刚度矩阵Stiffness(E,A,L)：输入弹性模量E，截面积A，长度L；输出单元刚度矩阵k
• 单元组装Assembly(KK,k,i,j)：输入单元刚度矩阵k，单元节点编号i、j；输出整体刚度矩阵KK
• 单元应力Stress(k,u,A)：输入单元矩阵k，单元位移矩阵u，截面积A；输出单元应力stress
• 单元节点力Force(k,u)：输入刚度矩阵k，单元位移矩阵u；输出节点力forces

2节点

1. 单元刚度矩阵函数

function k= Bar1D2Node_Stiffness(E,A,L)

k=[E*A/L -E*A/L;-E*A/L E*A/L];

2. 单元组装函数

function z= Bar1D2Node_Assembly(KK,k,i,j)

DOF(1)=i;
DOF(2)=j:
for nl=1:2
	for n2=1:2
		KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
	end
end

z=KK;

3. 单元应力函数

function stress= Bar1D2Node_Stress(k,u,A)

stress=k*E/A;

4. 单元节点力函数

function forces= Bar1D2Node_Force(k,u)

forces=k*u;

2D杆单元分析的matlab程序

设计四个函数

• 刚度矩阵Stiffness(E,A,x1,y1,x2,y2,alpha)：输入弹性模量E，截面积A，第一个节点坐标（x1,y1）,第二个节点坐标（x2,y2）和角度alpha（单位度）,；输出单元刚度矩阵k
• 单元组装Assembly(KK,k,i,j)：输入单元刚度矩阵k，单元节点编号i、j；输出整体刚度矩阵KK
• 单元应力Stress(E,x1,y1,x2,y2,alpha,u)：输入弹性模量E，第一个节点坐标（x1,y1）,第二个节点坐标（x2,y2），角度alpha（单位度），单元节点位移矩阵u；输出单元应力stress
• 单元节点力Force(E,A,x1,y1,x2,y2,alpha,u)：输入弹性模量E，截面积A，第一个节点坐标（x1,y1）,第二个节点坐标（x2,y2），角度alpha（单位度）和单元位移矩阵u；输出节点力forces

2个节点

1. 单元刚度矩阵函数

function k= Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha)

L=sqrt((x2-x1)*(x2-x1)+(y2-y1)(y2-y1))
x=alpha*pi/180
C=cos(x)
S=sin(x)
%单元刚度矩阵4*4
k=E*A/L[ C*C C*S -C*C -C*S;C*S S*S -C*S -S*S;
-C*C -C*S C*C C*S;-C*S -S*S C*S S*S];

2. 单元组装函数

function z= Bar1D2Node_Assembly(KK,k,i,j)

DOF(1)=2*i-1;
DOF(2)=2*i;
DOF(1)=2*j-1;
DOF(2)=2*j;
for nl=1:4
	for n2=1:4
		KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
	end
end

z=KK;

3. 单元应力函数

function stress= Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u)

L=sqrt((x2-x1)*(x2-x1)+(y2-y1)(y2-y1))
x=alpha*pi/180
C=cos(x)
S=sin(x)
stress=E/L*[-C -S C S]*u;

4. 单元节点力函数

function forces= Bar1D2Node_Force(k,u)

L=sqrt((x2-x1)*(x2-x1)+(y2-y1)(y2-y1))
x=alpha*pi/180
C=cos(x)
S=sin(x)
forces=E/L*[-C -S C S]*u*A;

四杆桁架的有限元分析的matlab编程

1. 结构的离散化与编号，见2.4章节

   1. E=2.95e11;A=0.0001;
   2. x1=0;y1=0;x2=0.4;y2=0;x3=0.4;y3=0.3;x4=0;y4=0.3;
   3. alphal=0;alpha2=90;alpha3=atan(0.75)*180/pi;

2. 计算各个单元的刚度矩阵

   • 分别针对4个杆单元，调用4次function k= Bar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha)，得到各个单元的刚度矩阵

3. 建立整体刚度矩阵方程

   • 该结构共有4个节点，设置整体的刚度矩阵为KK(8*8)
   • 先对KK清零，4次调用Bar1D2Node_Assembly(KK,k,i,j)
     

4. 边界条件的处理及刚度方程求解

   1. \(BC(u): u_1=0,v_1=0,v_2=0,u_4=0,v_4=0\)
   2. 总节点位移：\(\mathbf q=\begin{bmatrix} u_1 & v_1 & u_2 & v_2 & u_3 & v_3 & u_4 & v_4\end{bmatrix}^T\)
   3. 总节点力，R为支反力（‌支反力是指在物体受到外力作用时，在约束中产生的并作用在被约束物体上的力）：
      \(\mathbf P=\mathbf R+\mathbf F=\begin{bmatrix} R_{x1} & R_{y1} & 2\times10^4 & R_{y2} & 0 & -2.5\times10^4 & R_{x4} & R_{y4}\end{bmatrix}^T\)

5. 支反力的计算

6. 各个单元的应力计算

局部坐标系中的平面纯弯梁单元构建及MATLAB编程

平面纯弯梁的单元为2节点，每个节点具有2个自由度（位移和转角）

1. 节点描述

   • \(Node 1: x_1=0,Node 2: x_2=l\)
   • 节点位移矩阵：\(\mathbf q^e=\begin{bmatrix} v_1 & \theta_1 & v_2 & \theta_2\end{bmatrix}^T\)
   • 节点力移矩阵：\(\mathbf P^e=\begin{bmatrix} P_{v_1} & M_1 & P_{v_2} & M_2\end{bmatrix}^T\)
   • 刚度矩阵：\(\mathbf K^e \cdot \mathbf q^e =\mathbf P^e\)

2. 单元上的场描述——三大类变量

   • 位移场函数（挠度）- 使用多项式拟合：\(v(x)=a_0+a_1x+a_2x^2+a_3x^3\)
     节点条件：\(v|_{x=0}=v_1,v'|_{x=0}=\theta_1,v|_{x=l}=v_2,v'|_{x=0}=\theta_2\)
     解得，\(a_0=v_1,a_1=\theta _1\)，\(a_2=\frac{1}{l^2}(-3v_1-2\theta _1l+3v_2-\theta _2l)\)，\(a_3=\frac{1}{l^3}(2v_1+\theta _1l-2v_2+\theta _2l)\)

   • 令\(\xi _1 =x/l\)，则\(v(x)=(1-3\xi ^2+2\xi ^3)v_1+l(\xi -2\xi ^2+\xi ^3)\theta _1+(3\xi ^2-2 \xi ^3)v_2+l(\xi ^3- \xi ^2)\theta _2)\)

   • 最后的位移场函数：\(v(x)=\mathbf N(\xi) \cdot \mathbf q ^e\)

   • 形状函数矩阵：\(\mathbf N(\xi)=\begin{bmatrix} (1-3\xi ^2+2\xi ^3) & l(\xi -2\xi ^2+\xi ^3) & (3\xi ^2-2 \xi ^3) & l(\xi ^3- \xi ^2)\end{bmatrix}^T\)

   • 应变场函数：对挠度\(v(x)\)进行二阶求导：
     \(\epsilon (x,\hat{y})=-\hat{y} \frac{d^2v(x)}{dx^2}=\mathbf B(\xi)\cdot q^e\)，\(\hat y\)为相对于中心层为起点的y方向上的坐标

   • 几何矩阵：\(\mathbf B(\xi)=-\hat{y} \begin{bmatrix} B_1 & B_2 & B_3 & B_4\end{bmatrix}^T\)，其中\(B_1=\frac{1}{l^2}(12\xi -6),B_2=\frac{1}{l}(6\xi -4),B_3=-\frac{1}{l^2}(12\xi -6),B_4=\frac{1}{l}(6\xi -2)\)

   • 应力场方程：\(\sigma (x,\hat{y})=E \cdot \epsilon (x,\hat{y})=\mathbf S(\xi)\cdot \mathbf q^e\)，其中\(\mathbf S(\xi)=E^e \cdot \mathbf B(\xi)\)

3. 最小势能原理

   • 单元应变能：\(U^e=\frac{1}{2}\mathbf q^{eT}\cdot \mathbf K^e \cdot \mathbf q^e\)

   • 其中\(\mathbf K^e=\int _0^l \int _A \mathbf B^T \cdot E \cdot \mathbf B \cdot dA dx\)，\(\mathbf B(\xi)=-\hat{y} \begin{bmatrix} B_1 & B_2 & B_3 & B_4\end{bmatrix}^T\)

   • 单元的外力功：\(W^e=P_{v1}\cdot v_1+M_1\theta _1+P_{v2}\cdot v_2+M_2\theta _2=\mathbf P^{eT} \cdot \mathbf q^e\)

   • 应用最小势能原理
     

     

梁单元分析的matlab程序实现

5个函数

• 刚度矩阵Beam1D2Node_Stiffness(E,I,L)：该函数计算单元的刚度矩阵，输入弹性模量E，横截面的惯性矩I，梁单元的长度L，输出单元刚度矩阵k(4×4)
• 单元组装Beam1D2Node_Assembly(KK,k,i,j)：该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、i，输出整体刚度矩阵KK
• 单元应变Beam1D2Node_Strain(x,L,y)：该函数计算单元的几何矩阵，输入所测点距梁单元左节点的水平距离×，输入所测点以中性层为起点的y 方向的坐标，梁单元的长度L，输出单元几何形状函数矩阵B(1×4）
• 单元应力Beam1D2Node_Stress(E,B,u)：该函数计算单元内某点的应力，输入弹性模量E，几何矩阵B，节点位移列阵u，输出单元的应力stress
• 挠度BeamiD2NodeDeflection(x,Lu)：该函数计算单元内某点的挠度，输入所测点距梁单元左节点的水平距离x，梁单元的长度L，节点位移列阵u，输出该点的挠度v

function k=Beam1D2NodeStiffness(E,IL）
%该函数计算单元的刚度矩阵
%输入弹性模量E，横截面的惯性矩，梁单元的长度！
%输出单元刚度矩阵k(4×4）%

k = E*I/(L*L*L)*[12 6*L -12 6*L; 6*L 4*L*L -6*L 2*L*L;-12 -6*L 12 -6*L; 6*L 2*L*L -6*L 4*L*L];

function z=Beam1D2Node_Assembly(KK,ki,j)
%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k，单元的节点编号i、j
%输出整体刚度矩阵KK
%----------------------------------------

DOF(1)=2*i-1; 
DOF(2)=2*i; 
DOF(3)=2*j-1; 
DOF(4)=2*j; 
for n1=1:4 
	for n2=1:4
		KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2)
	end 
end 
z=KK;

function B=Beam1D2Node_Strain(x,L,y)
%该函数计算单元的几何矩阵
%输入所测点距梁单元左节点的水平距离×
%输入所测点以中性层为起点的y方向的坐标，梁单元的长度L
%输出单元几何形状函数矩阵B（1×4）%

e=x/L;
B1=(12*e-6)/(L*L); 
B2=(6*e-4)/L;
B3=-(12*e-6)/(L*L); 
B4=(6*e-2)/L;
B=-y*[B1,B2,B3,B4];

function stress=Beam1D2NodeStress(E,B,u）
%该函数计算单元内某点的应力
%输入弹性模量E，几何矩阵B，节点位移列阵u
%输出单元的应力stress

stress = E*B*u,

function v=Beam1D2Node_Deflection(x,L,u)
%该函数计算单元内某点的挠度
%输入所测点距梁单元左节点的水平距离x%输入梁单元的长度L，节点位移列阵u
%输出该点的挠度V

e=x/L;
Nl=1-3*e*e+2*e*e*e; 
N2=L(e-2*e*e+e*e*e); 
N3=3*e*e-2*e*e*e;
N4=L(e*e*e-e*e); 
N=[N1,N2,N3,N4]; 
v=N*u;

局部坐标系中的一般梁单元构建（组装）

1. 一般平面梁单元=平面纯弯梁单元+平面杆单元

2个节点单元

• 节点位移矩阵：\(\mathbf q^e=\begin{bmatrix} u_1 & v_1 & \theta_1 & u_2 & v_2 & \theta_2\end{bmatrix}^T\)
• 节点力移矩阵：\(\mathbf P^e=\begin{bmatrix} P_{u_1} & P_{v_1} & M_1 & P_{u_2} & P_{v_2} & M_2\end{bmatrix}^T\)
• 刚度矩阵：\(\mathbf K^e \cdot \mathbf q^e =\mathbf P^e\)

2. 一般空间梁单元

2个节点单元

• 节点位移矩阵：\(\mathbf q^e=\begin{bmatrix} u_1 & v_1 & w_1 & \theta_{x_1} & \theta_{y_1}& \theta_{z_1} & u_2 & v_2 & w_2 & \theta_{x_2}& \theta_{y_2}& \theta_{z_2}\end{bmatrix}^T\)
• 节点力移矩阵：\(\mathbf P^e=\begin{bmatrix} P_{u_1} & P_{v_1} & P_{w_1} & M_{x_1} & M_{y_1} & M_{z_1} & P_{u_2} & P_{v_2} & P_{w_2} & M_{x_2} & M_{y_2} & M_{z_2}\end{bmatrix}^T\)
• 刚度矩阵：\(\mathbf K^e \cdot \mathbf q^e =\mathbf P^e\)
• 其中\(\theta _x,\theta _y,\theta _z\)表示沿着全局坐标系相应轴旋转的角度

对应节点位移上的组装，得到12*12的刚度矩阵\(\mathbf K^e\)

梁单元的坐标变化

1. 平面梁单元的坐标变换

• 局部坐标系(ox)中的节点位移矩阵：\(\mathbf q^e=\begin{bmatrix} u_1 & v_1 & \theta_1 & u_2 & v_2 & \theta_2\end{bmatrix}^T\)
• 全局坐标系(\(\bar xo\bar y\) )中的节点力移矩阵：\(\bar {\mathbf q}^e=\begin{bmatrix} \bar u_1 & \bar v_1 & \theta_1 & \bar u_2 & \bar v_2 & \theta_2\end{bmatrix}^T\)
• 坐标转换方程：\(\mathbf q^e = \mathbf T^e \cdot \bar {\mathbf q}^e\)
• 局部坐标系刚度矩阵：\(\mathbf K^e \cdot \mathbf q^e =\mathbf P^e\)
• 全局坐标系刚度矩阵：\(\bar{\mathbf K}^e \cdot \bar{\mathbf q}^e =\bar{\mathbf P}^e\)
• 其中\(\bar{\mathbf K}^e={\mathbf T}^{eT} \cdot \mathbf K^e \cdot {\mathbf T}^{e}\)，\(\bar{\mathbf P}^e={\mathbf T}^{eT} \cdot \mathbf P^e\)， \({\mathbf T}^{e}\)如下图

2. 空间梁单元的坐标变换

• 局部坐标系(xoy)中的节点位移矩阵：
  \(\mathbf q^e=\begin{bmatrix} u_1 & v_1 & w_1 & \theta_{x_1} & \theta_{y_1}& \theta_{z_1} & u_2 & v_2 & w_2 & \theta_{x_2}& \theta_{y_2}& \theta_{z_2}\end{bmatrix}^T\)

• 全局坐标系(\(\bar xo\bar y\) )中的节点力移矩阵：
  \(\bar {\mathbf q}^e=\begin{bmatrix} \bar u_1 & \bar v_1 & \bar w_1 & \bar \theta_{x_1} & \bar \theta_{y_1}& \bar \theta_{z_1} & \bar u_2 & \bar v_2 & \bar w_2 & \bar \theta_{x_2}& \bar \theta_{y_2}& \bar \theta_{z_2}\end{bmatrix}^T\)

• 坐标系间的节点位移转换方程：\(\mathbf q^e = \mathbf T^e \cdot \bar {\mathbf q}^e\)

• 局部坐标系刚度矩阵：\(\mathbf K^e \cdot \mathbf q^e =\mathbf P^e\)

• 全局坐标系刚度矩阵：\(\bar{\mathbf K}^e \cdot \bar{\mathbf q}^e =\bar{\mathbf P}^e\)

• 其中\(\bar{\mathbf K}^e={\mathbf T}^{eT} \cdot \mathbf K^e \cdot {\mathbf T}^{e}\)，\(\bar{\mathbf P}^e={\mathbf T}^{eT} \cdot \mathbf P^e\)， \({\mathbf T}^{e}\)如下图

  

  

  

分布力的处理

• 集中力：作用在节点上的力，或者一个小面积上的力，该面积远小于任何方向上的尺寸，也可视为集中力
  • 支点力
  • 弯矩
• 分布力：当力的作用范围（面积）较大而不能忽略的作用力；有均布力和非均布力
  • 线分布力(\(N/m\))：楼板对梁的作用；符号q
  • 面分布力(\(N/m^2\))：风荷载对建筑物墙体的作用；符号p
  • 体分布力(\(N/m^3\))：物体的自重；符号v

使用外力功原理等效（不能使用静力学平衡等效），将单元上的分布力等效到结点的集中力

门型框架结构的实例分析及MATLAB编程

结构的参数：

• 弹性模量\(E=3.0\times 10^{11}Pa\)
• 惯性矩\(I = 6.5\times 10^{-7} m^4\)
• 杆截面积\(A = 6.8\times10^{-4}m^2\)
• 高\(H=0.96m\)

1. 结构离散化与编号

单元编号	节点	节点
1	1	2
2	3	1
3	4	2

按节点编号顺序

• 节点位移矩阵：
  \(\mathbf q=\begin{bmatrix} u_1 & v_1 & \theta_{1} &u_2 & v_2 & \theta_{2} &u_3 & v_3 & \theta_{3} &u_4 & v_4 & \theta_{4} \end{bmatrix}^T\)
• 节点外载矩阵：
  \(\mathbf F=\begin{bmatrix} F_{x_1} & F_{y_1} & M_{\theta_{1}}& 0 &F_{y_2} & M_{\theta _2} & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}^T\)
• 支反力矩阵：
  \(\mathbf R=\begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & R_{x_3} & R_{y_3} & R_{\theta_{3}}& R_{x_4} &R_{y_4} & R_{\theta_{4}} \end{bmatrix}^T\)
• 故总的节点载荷矩阵：
  \(\mathbf P=\mathbf F+\mathbf R=\)
  \(\begin{bmatrix} 3000 & -3000 & -720 & 0 & -3000 & 720 & R_{x_3} & R_{y_3} & R_{\theta_{3}}& R_{x_4} &R_{y_4} & R_{\theta_{4}} \end{bmatrix}^T\)

2. 各个单元的描述

单元1的局部坐标与整体坐标一致，故其刚度矩阵为：

单元2和单元3的情况相同，仅需要将节点编号进行映射即可

局部坐标系下的单元刚度矩阵为：

坐标转换矩阵为：

3. 组装整体刚度矩阵

\(\mathbf K=\mathbf K^{(1)}+\mathbf K^{(2)}+\mathbf K^{(3)}\)

\(\mathbf q= \sum \mathbf q^{(e)}\)，\(\mathbf P= \sum \mathbf P^{(e)}\)

刚度方程\(\mathbf K \cdot \mathbf q = \mathbf P\)

4. 边界条件

节点3、4固定：\(u_3=v_3=\theta _3 = u_4 = v_4 = \theta _4=0\)

处理边界条件后的刚度方程，只需要处理节点1和节点2

5. matlab编程

• 结构离散化与编号
• 各个单元的描述
  • 首先在MATLAB环境下，输入弹性模量E、横截面积A、惯性矩、长度L，由于单元2和单元3的刚度矩阵相同，针对单元1、单元2和单元3，分别二次调用函数 Beam2D2Node_Stiffness，就可以得到单元的刚度矩阵k1(6×6)和k2(6x6)。

>> E=3E11; 
>> I=6.5E-7; 
>> A=6.8E-4; 
>> L1=1.44; 
>> L2=0.96;
>> k1=Beam2D2Node_Stiffness(E,I,A,L1);
>> k2=Beam2D2Node_Stiffness(E,I,A,L2);

• 建立整体刚度矩阵：结构总的刚度矩阵KK为\(12*12\)，对KK置零，然后3次调用函数Beam2D2Node_Assemble进行刚度矩阵的组装

>> T=[0,1,0,0,0,0;-1,0,0,0,0,0;0,0,1,0,0,0;0,0,0,0,1,0;0,0,0,-1,0,0;0,0,0,0,0,1];
>> k3=T'*k2*T;
>> KK=zeros(12,12);
>> KK=Beam2D2Node_Assemble(KK,k1,1,2);
>> KK=Beam2D2Node_Assemble(KK,k3,3,1); 
>> KK=Beam2D2NodeAssemble(KK.k3.4.2);

• 边界条件

节点3、4固定：\(u_3=v_3=\theta _3 = u_4 = v_4 = \theta _4=0\)

>> k=KK(1:6,1:6);
>> p=[3000;-3000;-720;0;-3000;720];
>> u=k\p

先将上面得到的位移结果与位移边界条件的节点位移进行组合（注意位置关系），可以得到整体
的位移列阵U(12×1），再代回原整体刚度方程，计算出所有的节点力P(12×1）

>> U=[u;0;0;0;0;0;0]
U=0.0009 -0.0000 -0.0014 0.009 0 -0.000 -0.000
	0  0  0  0  0  0
>>P=KK*U
P=1.0e+003*
3.0000 -3.0000 -0.7200 0.0000 -3.0000 0.7200
-0.6658 2.2012 0.6014 -2.3342 3.7988 1.1283

门型框架结构的ANSYS实例分析

APDL，略

第8讲 连续体结构的有限元分析(1)

平面3节点三角形单元及MATLAB编程

平面三节点三角形单元，每个节点具有2个自由度，整体坐标系建模

1. 节点描述

• 节点坐标：\((x_i,y_i),i=1,2,3\)
• 节点位移矩阵：\(\mathbf q^e=\begin{bmatrix} u_1 & v_1 & u_2 & v_2 & u_3 & v_3\end{bmatrix}^T\)
• 节点力矩阵： \(\mathbf P ^e=\begin{bmatrix} P_{x_1} & P_{y_1} & P_{x_{2}}& P_{y_2} & P_{x_3} & P_{y_{3}} \end{bmatrix}^T\)

2. 单元上的场描述

• 位移场函数：所有节点间使用相同函数插值，低阶到高阶，唯一确定
  \(u_i(x,y)=\bar a_0+\bar a_1x_i+\bar a_2y_i\)
  \(v_i(x,y)=\bar b_0+\bar b_1x_i+\bar b_2y_i\)
  由\((u_i,v_i)\)可以确定系数\((\bar a_i,\bar b_i)\)的值，其中\(i=1,2,3\)
  \(u(x,y)=\bar a_0+\bar a_1x+\bar a_2y=N_1(x,y)u_1+N_2(x,y)u_2+N_3(x,y)u_3\)
  \(v(x,y)=\bar b_0+\bar b_1x+\bar b_2y=N_1(x,y)v_1+N_3(x,y)v_2+N_3(x,y)v_3\)
  \(N_i(x,y)=\frac{1}{2A}(a_i+b_ix+c_iy)\)
  例如\(a_1=x_2y_3-x_3y_2,b_1=y_2-y_3,c_1=-x_2+x_3,i=1,2,3\)

  位移场矩阵形式：\(\mathbf u(x,y)=\mathbf N(x,y) \cdot \mathbf q ^e\)
  其中\(\mathbf N(x,y)=\begin{bmatrix} N_1 & 0 & N_2 & 0 & N_3 & 0\\ 0 & N_1 & 0 & N_2 & 0 & N_3\end{bmatrix}\)

3. 应变场函数

\(\mathbf{\epsilon}(x,y)=\begin{bmatrix}\partial \end{bmatrix}_{3 \times 2} \mathbf{u}(x,y)=\begin{bmatrix}\partial \end{bmatrix}_{3 \times 2} \mathbf N(x,y) \cdot \mathbf q ^e=\mathbf B_{3\times 6}(x,y) \cdot \mathbf q ^e\)
其中，\(\begin{bmatrix}\partial \end{bmatrix}=\begin{bmatrix}\frac{\partial}{\partial x} & 0\\ 0 & \frac{\partial}{\partial y} \\\frac{\partial}{\partial y} & \frac{\partial}{\partial x} \end{bmatrix}\)，\(\mathbf B_{3\times 6}(x,y)=\begin{bmatrix}\partial \end{bmatrix}_{3 \times 2} \mathbf N(x,y)\)称为几何矩阵
令\(\mathbf B_{3\times 6}(x,y) =\begin{bmatrix} \mathbf (B_1)_{3\times 2}& \mathbf (B_2)_{3\times 2} & \mathbf (B_3)_{3\times 2}\end{bmatrix}\)
那么，\(\mathbf (B_i)_{3\times 2}=\frac{1}{2A} \begin{bmatrix} b_i & 0 \\ 0 & c_i \\ c_i & b_i \end{bmatrix},i=1,2,3\)

• 应力场函数
  \(\mathbf {\sigma}(x,y,z)_{(3\times 1)}=\frac{E}{1-\mu ^2} \begin{bmatrix} 1 & \mu & 0\\ \mu & 1 & 0\\ 0& 0 & \frac{1-\mu}{2} \end{bmatrix}\begin{bmatrix}\epsilon _{xx} \\\epsilon _{yy} \\ \gamma _{xy}\end{bmatrix}=\mathbf{ D_{(3\times 3)}} \cdot \epsilon=\mathbf{ D} \cdot (\mathbf B \cdot \mathbf q ^e)\)
  其中，\(\mathbf D\)称为弹性系数矩阵
  令\(\mathbf S=\mathbf D \cdot \mathbf B\)，并称\(\mathbf S\)为应力函数矩阵

4. 最小势能原理
   \(\prod ^e=U-W=\frac{1}{2} \int _{\Omega ^e} \sigma _{ij} \epsilon _{ij} d \Omega-[\int _{\Omega ^e} \bar b_i u _i d \Omega+\int _{S_p ^e} \bar p_i u _i dA]\)
   \(=\frac{1}{2} \mathbf q^{eT}(\int _{\Omega ^e} \mathbf B^{T}\mathbf D \mathbf B d \Omega)\mathbf q^{e}-(\int _{\Omega ^e} \mathbf N^{T} \mathbf{\bar b} d \Omega+\int _{S_p ^e} \mathbf N^{T} \mathbf{\bar p} dA)^T \mathbf q^{e}\)
   \(=\frac{1}{2} \mathbf q^{eT} \mathbf K ^e\mathbf q^{e}-\mathbf P^{eT} \mathbf q^{e}\)
   其中，\(\mathbf K ^e_{(6\times 6)}=\int _{\Omega ^e} \mathbf B^{T}\mathbf D_{(3\times 3)} \mathbf B d \Omega=\int _{\Omega ^e} \mathbf B^{T}\mathbf D \mathbf B dA \cdot t\)，每个元素为
   
   \(\mathbf P^{eT} =\int _{\Omega ^e} \mathbf N^{T} \mathbf{\bar b} d \Omega+\int _{S_p ^e} \mathbf N^{T} \mathbf{\bar p} dA\)

应用最小势能原理，\(\frac{\partial \prod(\mathbf q^{e}) }{\partial\mathbf q^{e}}=0\)
则获得刚度方程，\(\mathbf K ^e \cdot \mathbf q^{e}=\mathbf P^{e}\)

• 平面三节点三角形单元性质
  • 单元内的应力和应变均为常数，又称为常应变三角形单元CST(constant strain triangle)
  • 没有针对节点位移的坐标变换
  • 对于应变梯度较大的区域，CST单元划分应适当加密，否则导致较大误差

5. 平面三节点三角形单元的matlab编程

   三个函数

• 刚度矩阵Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj.yj,xm,ym,ID)：该函数计算单元的刚度矩阵，输入弹性模量E，泊松比NU，厚度t，三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym，平面问题性质指示参数ID(1为平面应力，2为平面应变），输出单元刚度矩阵k(6X6)
• 单元组装Triangle2D3Node_Assembly(Kk,k,i,j,m)：该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、j、m，输出整体刚度矩阵KK
• 单元应力Triangle2D3Node_Stress(E,NU,xi,yi,xj.yj,xm,ym,u,ID）：该函数计算单元的应力，输入弹性模量E，泊松比NU，厚度t，三个节点i、j、m的坐标xi,yi,xj,yj,xm,ym，平面问题性质指标参数ID(1为平面应力，2为平面应变），单元的位移列阵u(6X1），输出单元的应力stress，由于它为常应力单元，则单元的应力分量为Sx,Sy,Sxy

function k=Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)
%该函数计算单元的刚度矩阵
%输入弹性模量E，泊松比NU，厚度t
%输入三个节点i、j、m的坐标xi,yi,xj,yi,xm,ym
%输入平面问题性质指标参数ID(1为平面应力，2为平面应变）
%输出单元刚度矩阵k(6X6)

A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yi))/2;
betai = yj-ym; 
betaj = ym-yi; 
betam = yi-yj; 
gammai=xm-xj; 
gammaj=xi-xm; 
gammam=xj-xi;
B=[betai 0 betaj 0 betam 0;
0 gammai 0 gammaj 0 gammam;
gammai betai gammaj betaj gammam betam]/(2*A); 
if ID == 1
	D=(E/(1-NU*NU)*[1 NU 0;NU 1 0;0 0 (1-NU)/2];
elseif ID == 2
	D=(E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2];
end
k=t*A*B'*D*B;

function z=Triangle2D3Node_Assembly(KK,k,i,j,m)
%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k
%输入单元的节点编号i、j、m
%输出整体刚度矩阵KK

DOF(1)=2*i-1; 
DOF(2)=2*i; 
DOF(3)=2*j-1; 
DOF(4)=2*j; 
DOF(5)=2*m-1; 
DOF(6)=2*m; 
for n1=1:6 
	for n2=1:6
		KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2);
	end 
end 
z=KK;

function stress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)
%该函数计算单元的应力
%输入弹性模量E，泊松比NU，厚度t
%输入三个节点i、j、m的坐标xi,yi,xj.yj,xm,ym
%输入平面问题性质指标参数ID(1为平面应力，2为平面应变），单元的位移列阵u(6X1）
%输出单元的应力stress(3X1），由于它为常应力单元，则单元的应力分量为Sx,Sy,Sxy

A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yi))/2;
betai = yj-ym; 
betaj = ym-yi; 
betam=yi-yi 
gammai=xm-xj; 
gammaj=xi-xm;
gammam=xj-xi;
B=[betai 0 betaj 0 betam 0;
0 gammai 0 gammaj 0 gammam;
gammai betai gammaj betaj gammam betam]/(2*A); 
if ID == 1
	D = (E/(1-NU*NU)*[1 NU 0; NU 1 0;0 0 (1-NU)/2]; 
elseif ID==2
	D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2]; 
end
stress = D*B*u;

平面4节点矩形单元及MATLAB编程

平面4节点矩形单元，每个节点具有2个自由度，局部坐标系建模

令\(\xi=\frac{x}{a},\eta =\frac{y}{b}\)

1. 节点描述

• 节点坐标(逆时针编号)：\((x_i,y_i),i=1,2,3,4\)
• 节点位移矩阵：\(\mathbf q^e=\begin{bmatrix} u_1 & v_1 & u_2 & v_2 & u_3 & v_3 & u_4 & v_4\end{bmatrix}^T\)
• 节点力矩阵： \(\mathbf P ^e=\begin{bmatrix} P_{x_1} & P_{y_1} & P_{x_{2}}& P_{y_2} & P_{x_3} & P_{y_{3}}& P_{x_4} & P_{y_{4}} \end{bmatrix}^T\)

2. 位移场函数

• 位移场函数：所有节点间使用相同函数插值，低阶到高阶，唯一确定
  \(u_i(x,y)= a_0+ a_1x_i+ a_2y_i+a_3x_iy_i\)
  \(v_i(x,y)= b_0+ b_1x_i+ b_2y_i+b_3x_iy_i\)
  由\((u_i,v_i)\)可以确定系数\(( a_i, b_i)\)的值，其中\(i=1,2,3,4\)

  \[u(x,y)= a_0+ a_1x+ a_2y+a_3xy=N_1(x,y)u_1+N_2(x,y)u_2+N_3(x,y)u_3+N_4(x,y)u_4 \]
  \[v(x,y)= b_0+ b_1x+ b_2y+b_3xy=N_1(x,y)v_1+N_3(x,y)v_2+N_3(x,y)v_3+N_4(x,y)v_4 \]

  无量纲\(N_i=\frac{1}{4}(1+\xi _i \xi)(1+\eta _i \eta)\)为形状函数
  例如
  

  位移场矩阵形式：\(\mathbf u(x,y)=\mathbf N(x,y) \cdot \mathbf q ^e\)
  其中\(\mathbf N(x,y)=\begin{bmatrix} N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 & 0\\ 0 & N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4\end{bmatrix}\)

3. 应变场函数
   \(\mathbf{\epsilon}(x,y)=\begin{bmatrix}\partial \end{bmatrix}_{3 \times 2} \mathbf{u}(x,y)=\begin{bmatrix}\partial \end{bmatrix}_{3 \times 2} \mathbf N(x,y) \cdot \mathbf q ^e=\mathbf B_{3\times 8}(x,y) \cdot \mathbf q ^e\)
   其中，\(\begin{bmatrix}\partial \end{bmatrix}=\begin{bmatrix}\frac{\partial}{\partial x} & 0\\ 0 & \frac{\partial}{\partial y} \\\frac{\partial}{\partial y} & \frac{\partial}{\partial x} \end{bmatrix}\)，\(\mathbf B_{3\times 8}(x,y)=\begin{bmatrix}\partial \end{bmatrix}_{3 \times 2} \mathbf N(x,y)\)称为几何矩阵
   令\(\mathbf B_{3\times 8}(x,y) =\begin{bmatrix} \mathbf (B_1)_{3\times 2}& \mathbf (B_2)_{3\times 2} & \mathbf (B_3)_{3\times 2}& \mathbf (B_4)_{3\times 2}\end{bmatrix}\)
   那么，\(\mathbf (B_i)_{3\times 2}= \begin{bmatrix} \frac{\partial N_i}{\partial x} & 0 \\ 0 & \frac{\partial N_i}{\partial y} \\ \frac{\partial N_i}{\partial y} & \frac{\partial N_i}{\partial x} \end{bmatrix},i=1,2,3,4\)

• 应力场函数
  \(\mathbf {\sigma}(x,y,z)_{(3\times 1)}=\frac{E}{1-\mu ^2} \begin{bmatrix} 1 & \mu & 0\\ \mu & 1 & 0\\ 0& 0 & \frac{1-\mu}{2} \end{bmatrix}\begin{bmatrix}\epsilon _{xx} \\\epsilon _{yy} \\ \gamma _{xy}\end{bmatrix}=\mathbf{ D_{(3\times 3)}} \cdot \epsilon=\mathbf{ D} \cdot (\mathbf B \cdot \mathbf q ^e)\)
  其中，\(\mathbf D\)称为弹性系数矩阵
  令\(\mathbf S=\mathbf D \cdot \mathbf B\)，并称\(\mathbf S\)为应力函数矩阵

4. 最小势能原理
   \(\prod ^e=U-W=\frac{1}{2} \int _{\Omega ^e} \sigma _{ij} \epsilon _{ij} d \Omega-[\int _{\Omega ^e} \bar b_i u _i d \Omega+\int _{S_p ^e} \bar p_i u _i dA]\)
   \(=\frac{1}{2} \mathbf q^{eT}(\int _{\Omega ^e} \mathbf B^{T}\mathbf D \mathbf B d \Omega)\mathbf q^{e}-(\int _{\Omega ^e} \mathbf N^{T} \mathbf{\bar b} d \Omega+\int _{S_p ^e} \mathbf N^{T} \mathbf{\bar p} dA)^T \mathbf q^{e}\)
   \(=\frac{1}{2} \mathbf q^{eT} \mathbf K ^e\mathbf q^{e}-\mathbf P^{eT} \mathbf q^{e}\)
   其中，单元刚度矩阵\(\mathbf K ^e_{(8\times 8)}=\int _{\Omega ^e} \mathbf B^{T}\mathbf D_{(3\times 3)} \mathbf B d \Omega=\int _{\Omega ^e} \mathbf B^{T}\mathbf D \mathbf B dA \cdot t\)，每个元素为

\(\mathbf P^{eT} =\int _{\Omega ^e} \mathbf N^{T} \mathbf{\bar b} d \Omega+\int _{S_p ^e} \mathbf N^{T} \mathbf{\bar p} dA\)

应用最小势能原理，\(\frac{\partial \prod(\mathbf q^{e}) }{\partial\mathbf q^{e}}=0\)
则获得刚度方程，\(\mathbf K ^e \cdot \mathbf q^{e}=\mathbf P^{e}\)

• 平面四节点矩形单元性质
  • 位移场函数=完全线性+交叉项函数，属于双线性位移模式
  • 应变场和应力场函数均为不完全线性函数
  • 平面4节点矩形单元比3节点三角形单元的精度高
  • 矩形单元的几何形状还应变换为任意边形的形状，以增强几何的适应性

5. 平面三节点三角形单元的matlab编程

   三个函数

• 刚度矩阵Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj.yj,xm,ym,xp,yp,ID)：
  该函数计算单元的刚度矩阵，输入弹性模量E，泊松比NU，厚度h，4个节点i、j、m的坐标xi,yi,xj,yj,xm,ym,xp,yp，平面问题性质指示参数ID(1为平面应力，2为平面应变），输出单元刚度矩阵k(8X8)
• 单元组装Quad2D4Node_Assembly(Kk,k,i,j,m,p)：该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、j、m、p，输出整体刚度矩阵KK
• 单元应力Quad2D4Node_Stress(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,u,ID）：
  该函数计算单元的应力，输入弹性模量E，泊松比NU，厚度h，4个节点i、j、m、p的坐标xi,yi,xj,yj,xm,ym,xp,yp，平面问题性质指标参数ID(1为平面应力，2为平面应变），单元的位移列阵u(8X1），输出单元的应力stress，由于它为常应力单元，则单元的应力分量为Sx,Sy,Sxy

function k=Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj.yj,xm,ym,xp.yp,ID)
%该函数计算单元的刚度矩阵
%输入弹性模量E，泊松比NU，厚度h
%输入4个节点i、j、m、p的坐标xi,yi,xj,yi,xm,ym,xp,yp
%输入平面问题性质指标参数ID(1为平面应力，2为平面应变）
%输出单元刚度矩阵k(8X8)

syms s t;
a = (yi*(s-1)+yi*(-1-s)+ym*(1+s)+yp*(1-s))/4; 
b = (yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4; 
c = (xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4; 
d = (xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s))/4;

B1=[a*(t-1)/4-b*(s-1)/4 0;0 c*(s-1)/4-d*(t-1)/4;c*(s-1)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4];
B2=[a*(1-t)/4-b*(-1-s)/4 0;0 c*(-1-s)/4-d*(1-t)/4;c*(-1-s)/4-d*(1-t)/4 a*(1-t)/4-b*(-1-s)/4];
B3=[a*(t+1)/4-b*(s+1)/4 0;0 c*(s+1)/4-d*(t+1)/4;c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4]
B4=[a*(-1-t)/4-b*(1-s)/4 0;0 c*(1-s)/4-d*(-1-t)/4;c*(1-s)/4-d*(-1-t)/4 a*(-1-t)/4-b*(1-s)/4];

Bfirst=[B1 B2 B3 B4];
Jfirst =[0 1-t t-s s-1;t-1 0 s+1 -s-t;s-t -s-1 0 t+1;1-s s+t -t-1 o];
J=[xi xj xm xp]*Jfirst*[yi;yj;ym;yp]/8;
B=Bfirst/J; 
if ID == 1
	D=(E/(1-NU*NU))*[1 NU 0;NU 1 0;0 0 (1-NU)/2];
elseif ID == 2
	D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2];
end

BD=J*transpose(B)*D*B;
r = int(int(BD, t, -1, 1), s, -1, 1); 
z = h*r;
k = double(z);

function z=Quad2D4Node_Assembly(KK,k,i,j.m,p)
%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k，单元的节点编号i、j、m、p
%输出整体刚度矩阵KK

DOF(1)=2*i-1; 
DOF(2)=2*i; 
DOF(3)=2*j-1; 
DOF(4)=2*j; 
DOF(5)=2*m-1; 
DOF(6)=2*m; 
DOF(7)=2*p-1; 
DOF(8)=2*p; 
for n1=1:8 
	for n2=1:8
		KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
	end 
end 
z=KK;

function stress=Quad2D4Node_Stress(E,NU,xi,yi,xj.yj,xm,ym,xp.yp,u,ID)
%该函数计算单元的应力
%输入弹性模量E，泊松比NU，厚度h
%输入4个节点i、j、m、p的坐标xi,yi,xj.yj,xm,ym,xp,yp，
%输入平面问题性质指标参数ID(1为平面应力，2为平面应变）
%输入单元的位移列阵u(8X1)
%输出单元的应力stress(3x1)
%由于它为常应力单元，则单元的应力分量为Sx,Sy,Sxy

syms st;
a = (yi*(s-1)+yi*(-1-s)+ym*(1+s)+yp*(1-s))/4; 
b = (yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t)/4; 
c = (xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t)/4;
d = (xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s)/4;

B1 = [a*(t-1)/4-b*(s-1)/4 0;0 c*(s-1)/4-d*(t-1)/4;c*(s-1)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4];
B2 =[a*(1-t)/4-b*(-1-s)/4 0;0 c*(-1-s)/4-d*(1-t)/4;c*(-1-s)/4-d*(1-t)/4 a*(1-t)/4-b*(-1-s)/4];
B3 = [a*(t+1)/4-b*(s+1)/4 0;0 c*(s+1)/4-d*(t+1)/4;c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4];
B4 = [a*(-1-t)/4-b*(1-s)/4 0;0 c*(1-s)/4-d*(-1-t)/4;c*(1-s)/4-d*(-1-t)/4 a*(-1-t)/4-b*(1-s)/4];

Bfirst = [B1 B2 B3 B4];
Jfirst = [0 1-t t-s s-1;t-1 0 s+1 -s-t;s-t -s-1 0 t+1;1-s s+t -t-1 0];
J=[xi xj xm xp]*Jfirst*[yi;yj;ym;yp]/8;
B=Bfirst/J; 
if ID == 1
	D = (E/(1-NU*NU))*[1 NU 0; NU 1 0;0 0 (1-NU)/2];
elseif ID == 2
	D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0;NU 1-NU 0;0 0 (1-2*NU)/2]
end
str1=D*B*u;
str2= subs(str1,{s,t}，{0,0}); 
stress=double(str2);

轴对称单元

轴对称物体的单元采用柱坐标表述，旋转截面离散成平面单元

三大类力学变量

• 位移：\(\begin{bmatrix} u_r & w\end{bmatrix}^T,u_\theta =0\)
• 应变：\(\begin{bmatrix} \epsilon _{rr} & \epsilon _{\theta \theta } & \epsilon _{zz} & \gamma _{rz}\end{bmatrix}^T,\gamma _{r\theta }=\gamma _{\theta z} =0\)
• 应力：\(\begin{bmatrix} \sigma _{rr} & \sigma _{\theta \theta } & \sigma _{zz} & \tau _{rz}\end{bmatrix}^T,\tau _{r\theta }=\tau _{\theta z} =0\)

三大类方程

使用三节点三角形单元

1. 节点描述

• 节点坐标：\((r_i,z_i),i=1,2,3\)
• 节点位移矩阵：\(\mathbf q^e=\begin{bmatrix} u_{r_1} & w_1 & u_{r_2} & w_2 & u_{r_3} & w_3\end{bmatrix}^T\)
• 节点力矩阵： \(\mathbf P ^e=\begin{bmatrix} P_{r_1} & P_{z_1} & P_{r_{2}}& P_{z_2} & P_{r_3} & P_{z_{3}} \end{bmatrix}^T\)

2. 单元上的场描述

• 位移场函数：所有节点间使用相同函数插值，低阶到高阶，唯一确定
  \(u_{r_i}(r,z)=\bar a_0+\bar a_1r_i+\bar a_2z_i\)
  \(w_i(r,z)=\bar b_0+\bar b_1r_i+\bar b_2z_i\)
  由\((u_{r_i},w_i)\)可以确定系数\((\bar a_i,\bar b_i)\)的值，其中\(i=1,2,3\)
  \(u(r,z)=\bar a_0+\bar a_1r+\bar a_2z=N_1(r,z)u_{r_1}+N_2(r,z)u_{r_2}+N_3(r,z)u_{r_3}\)
  \(w(r,z)=\bar b_0+\bar b_1r+\bar b_2z=N_1(r,z)w_1+N_3(r,z)w_2+N_3(r,z)w_{3}\)
  \(N_i(r,z)=\frac{1}{2A}(a_i+b_ir+c_iz)\)
  例如\(a_1=r_2z_3-r_3z_2,b_1=z_2-z_3,c_1=-r_2+r_3,i=1,2,3\)

  位移场矩阵形式：\(\mathbf u(r,z)=\mathbf N(r,z) \cdot \mathbf q ^e\)
  其中\(\mathbf N(r,z)=\begin{bmatrix} N_1 & 0 & N_2 & 0 & N_3 & 0\\ 0 & N_1 & 0 & N_2 & 0 & N_3\end{bmatrix}\)

3. 应变场函数

\(\mathbf{\epsilon}(r,z)=\begin{bmatrix}\partial \end{bmatrix}_{4 \times 2} \mathbf{u}(r,z)=\begin{bmatrix}\partial \end{bmatrix}_{4 \times 2} \mathbf N(r,z) \cdot \mathbf q ^e=\mathbf B_{4\times 6}(r,z) \cdot \mathbf q ^e\)
其中，\(\begin{bmatrix}\partial \end{bmatrix}=\begin{bmatrix}\frac{\partial}{\partial r} & 0\\ \frac{1}{r} & 0 \\ 0 & \frac{\partial}{\partial z} \\ \frac{\partial}{\partial z} & \frac{\partial}{\partial r} \end{bmatrix}\)，\(\mathbf B_{4\times 6}(r,z)=\begin{bmatrix}\partial \end{bmatrix}_{4 \times 2} \mathbf N(r,z)\)称为几何矩阵
令\(\mathbf B_{3\times 6}(r,z) =\begin{bmatrix} \mathbf (B_1)_{3\times 2}& \mathbf (B_2)_{3\times 2} & \mathbf (B_3)_{3\times 2}\end{bmatrix}\)
那么，\(\mathbf (B_i)_{3\times 2}=\frac{1}{2A} \begin{bmatrix} b_i & 0 \\ 0 & c_i \\ c_i & b_i \end{bmatrix},i=1,2,3\)

• 应力场函数
  \(\mathbf {\sigma}_{(4\times 1)}=\mathbf{ D_{(4\times 4)}} \cdot \epsilon=\mathbf{ D} \cdot (\mathbf B \cdot \mathbf q ^e)=\mathbf S \cdot \mathbf q ^e\)
  其中，\(\mathbf D\)称为弹性系数矩阵

  \[\mathbf D=\frac{E(1-\mu )}{(1+\mu)(1-2\mu )} \begin{bmatrix} 1 & \frac{\mu}{1-\mu} & \frac{\mu}{1-\mu} & 0\\ \frac{\mu}{1-\mu} & 1 & \frac{\mu}{1-\mu} & 0\\ \frac{\mu}{1-\mu} & \frac{\mu}{1-\mu} & 1 & 0\\ 0 & 0 & 0 & \frac{1-2\mu}{2(1-\mu)}\end{bmatrix} \]

  \(\mathbf S=\mathbf D \cdot \mathbf B\)，并称\(\mathbf S\)为应力函数矩阵

4. 最小势能原理

\(\prod ^e=U-W\)
\(=\frac{1}{2} \mathbf q^{eT} \mathbf K ^e\mathbf q^{e}-\mathbf P^{eT} \mathbf q^{e}\)
其中，
\(\mathbf K ^e_{(6\times 6)}=\int _{\Omega ^e} \mathbf B^{T}\mathbf D_{(4\times 4)} \mathbf B d \Omega=\int _{A ^e}\int _{0}^{2\pi} \mathbf B^{T}\mathbf D \mathbf B \cdot rdrd\theta dz=2\pi r\int _{A ^e}\mathbf B^{T}\mathbf D \mathbf B \cdot dr dz\)

\(\mathbf P^{eT} =\int _{\Omega ^e} \mathbf N^{T} \mathbf{\bar b} d \Omega+\int _{S_p ^e} \mathbf N^{T} \mathbf{\bar p} dA=2\pi r(\int _{A ^e} \mathbf N^{T} \mathbf{\bar b} drdz+\int _{l_p ^e} \mathbf N^{T} \mathbf{\bar p} dl)\)

应用最小势能原理，\(\frac{\partial \prod(\mathbf q^{e}) }{\partial\mathbf q^{e}}=0\)
则获得刚度方程，\(\mathbf K ^e \cdot \mathbf q^{e}=\mathbf P^{e}\)

四节点举行环形单元

1. 节点描述

• 节点坐标：\((r_i,z_i),i=1,2,3,4\)
• 节点位移矩阵：\(\mathbf q^e=\begin{bmatrix} u_{r_1} & w_1 & u_{r_2} & w_2 & u_{r_3} & w_3& u_{r_4} & w_4\end{bmatrix}^T\)
• 节点力矩阵： \(\mathbf P ^e=\begin{bmatrix} P_{r_1} & P_{z_1} & P_{r_{2}}& P_{z_2} & P_{r_3} & P_{z_{3}} & P_{r_4} & P_{z_{4}} \end{bmatrix}^T\)

2. 单元上的场描述

• 位移场函数：所有节点间使用相同函数插值，低阶到高阶，唯一确定
  \(u_{r_i}(r,z)= a_0+ a_1r_i+ a_2z_i+a_3r_iz_i\)
  \(w_i(r,z)= b_0+ b_1r_i+ b_2z_i+b_3r_iz_i\)

2. 最小势能原理

\(\prod ^e=U-W\)
\(=\frac{1}{2} \mathbf q^{eT} \mathbf K ^e\mathbf q^{e}-\mathbf P^{eT} \mathbf q^{e}\)

其中，
\(\mathbf K ^e_{(8\times 8)}=\int _{\Omega ^e} \mathbf B^{T}\mathbf D_{(4\times 4)} \mathbf B d \Omega=\int _{A ^e}\int _{0}^{2\pi} \mathbf B^{T}\mathbf D \mathbf B \cdot rdrd\theta dz=2\pi r\int _{A ^e}\mathbf B^{T}\mathbf D \mathbf B \cdot dr dz\)

\(\mathbf P^{eT} =\int _{\Omega ^e} \mathbf N^{T} \mathbf{\bar b} d \Omega+\int _{S_p ^e} \mathbf N^{T} \mathbf{\bar p} dA=2\pi r(\int _{A ^e} \mathbf N^{T} \mathbf{\bar b} drdz+\int _{l_p ^e} \mathbf N^{T} \mathbf{\bar p} dl)\)

应用最小势能原理，\(\frac{\partial \prod(\mathbf q^{e}) }{\partial\mathbf q^{e}}=0\)
则获得刚度方程，\(\mathbf K ^e \cdot \mathbf q^{e}=\mathbf P^{e}\)

分布力的处理

需要将分布力转换等效为节点载荷

最小势能原理

\(\prod ^e=U^e-W^e=\frac{1}{2} \int _{\Omega ^e} \sigma _{ij} \epsilon _{ij} d \Omega-W^e\)

其中，\(\mathbf u = \mathbf N \cdot \mathbf q^e\)

\(W^e=\int _{\Omega ^e} \bar b_i u _i d \Omega+\int _{S_p ^e} \bar p_i u _i dA=\int _{\Omega ^e} \mathbf{\bar b}^T \mathbf u d \Omega+\int _{S_p ^e} \mathbf{\bar p} ^T\mathbf u dA=\mathbf P^{eT} \cdot \mathbf q^e\)

等效节点载荷（平面三节点三角形单元）

平面矩形薄板分析的MATLAB编程

问题描述

离散成两个3节点三角形单元，外力载荷离散到节点力

结构的参数：

• 弹性模量\(E=10^{7}Pa\)
• 载荷\(F=100000N\)
• 泊松比\(\mu = 1/3\)
• 厚度\(t=0.1m\)

在MATLAB环境下，输入弹性模量E、泊松比NU，薄板厚度为t，平面应力问题性质指示参数ID，然后针对单元1和单元2，分别调用两次函数Triangle2D3Node_Stiffness，就可以得到单元的刚度矩阵k1(6×6)和k2(6×6)

窗口输入

>> E=1e7; 
>> NU=1/3;  
>> t=0.1;  
>> ID=1; 
>> k1=Triangle2D3Node_Stiffness(E,NU,t,2,0,0,1,0,0,ID)
>> k2=Triangle2D3Node_Stiffness(E,NU,t,0,1,2,0,2,1,ID)

由于该结构共有4个节点，则总共的自由度数为8，因此，结构总的刚度矩阵为KK(8×8)，先对KK清零，然后两次调用函数Triangle2D3Node_Assemblv进行刚度矩阵的组装

>> kk = zeros(8,8);
>> KK=Triangle2D3Node_Assembly(KK,k1,2,3,4) 
>> KK=Triangle2D3NodeAssembly(KK.k2.3.2.1)

边界条件，\(u_3=0,v_3=0,u_4=0,v_4=0\)

将针对节点1和节点2的位移进行求解，节点1和节点2的位移将对应KK矩阵中的前4行和前4列，则需从KK(8×8)中提出，置给k，然后生成对应的载荷列阵p，再采用高斯消去法进行求解。
注意：MATLAB中的"\"就采用用高斯消去法

>> k = kk(1:4,1:4);
>> p = [0;-50000;0;-50000] 
>> u=k\p
u = 0.188 -0.899 -0.150 -0.842

即，\(u_1=0.188,v_1=-0.899,u_2=-0.151,v_2=-0.842\)

在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力

>> U=[u;0;0;0;0]; 
>> P=KK*U
P=1.0e+005 *
-0.0000 -0.5000 0 -0.5000 -2.0000 -0.0702 2.0000 1.0702

即，\(R_{x_3}=-200000,R_{y_3}=-7020,R_{x_4}=200000,R_{y_4}=107020\)

平面矩形薄板的ANSYS实例分析

第9讲 连续体结构的有限元分析(2)

空间4节点四面体单元及MATLAB编程

每个节点具有3个自由度，整体坐标系建模

1. 节点描述

• 节点坐标：\((x_i,y_i,z_i),i=1,2,3,4\)
• 节点位移矩阵：
  \(\mathbf q^e=\begin{bmatrix} u_1 & v_1 & w_1 & u_2 & v_2 & w_2 & u_3 & v_3 & w_3 & u_4 & v_4 & w_4\end{bmatrix}^T\)
• 节点力矩阵：
  \(\mathbf P ^e=\begin{bmatrix} P_{x_1} & P_{y_1} & P_{z_1}& P_{x_{2}}& P_{y_2}& P_{z_2} & P_{x_3} & P_{y_{3}} & P_{z_3} & P_{x_4} & P_{y_{4}} & P_{z_4} \end{bmatrix}^T\)

2. 单元上的场描述

• 位移场函数：所有节点间使用相同函数插值，低阶到高阶，唯一确定
  \(u_i(x,y,z)=\bar a_0+\bar a_1x_i+\bar a_2y_i+\bar a_3z_i\)
  \(v_i(x,y,z)=\bar b_0+\bar b_1x_i+\bar b_2y_i+\bar b_3z_i\)
  \(w_i(x,y,z)=\bar c_0+\bar c_1x_i+\bar c_2y_i+\bar c_3z_i\)
  由\((u_i,v_i,w_i)\)可以确定系数\((\bar a_i,\bar b_i,\bar c_i)\)的值，其中\(i=1,2,3,4\)
  \(u(x,y,z)=\bar a_0+\bar a_1x+\bar a_2y+\bar a_3z=N_1(x,y,z)u_1+N_2u_2+N_3u_3+N_4u_4\)
  \(v(x,y,z)=\bar b_0+\bar b_1x+\bar b_2y+\bar b_3z=N_1(x,y,z)v_1+N_3v_2+N_3v_3+N_3v_4\)
  \(w(x,y,z)=\bar c_0+\bar c_1x+\bar c_2y+\bar c_3z=N_1(x,y,z)w_1+N_2w_2+N_3w_3+N_4w_4\)

\(N_i(x,y,z)=\frac{1}{6V}(a_i+b_ix+c_iy+d_iz),i=1,2,3,4\)

位移场矩阵形式：\(\mathbf u(x,y,z)=\mathbf N(x,y,z) \cdot \mathbf q ^e\)
其中\(\mathbf N(x,y,z)=\begin{bmatrix} N_1 & 0 & 0 & N_2 & 0& 0 & N_3 & 0 & 0 & N_4 & 0 & 0\\ 0 & N_1 & 0 & 0 & N_2 & 0 & 0 & N_3 & 0 & 0 & N_4& 0 \\ 0 & 0 & N_1 & 0 & 0 &N_2 &0 &0 &N_3& 0 &0 &N_4 \end{bmatrix}\)

3. 应变场函数

\(\mathbf{\epsilon}(x,y,z)=\begin{bmatrix}\partial \end{bmatrix}_{6 \times 3} \mathbf{u}(x,y)=\begin{bmatrix}\partial \end{bmatrix}_{6 \times 3} \mathbf N(x,y) \cdot \mathbf q ^e=\mathbf B_{6\times 12}(x,y) \cdot \mathbf q ^e\)
其中，\(\begin{bmatrix}\partial\end{bmatrix}=\begin{bmatrix} \frac{\partial}{\partial x} & 0 & 0\\ 0 & \frac{\partial}{\partial y} & 0\\0 & 0 & \frac{\partial}{\partial z}\\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0\\ 0 & \frac{\partial}{\partial z}& \frac{\partial}{\partial y}\\ \frac{\partial}{\partial z} & 0 & \frac{\partial}{\partial x}\end{bmatrix}\)，\(\mathbf B_{6\times 12}(x,y,z)=\begin{bmatrix}\partial \end{bmatrix}_{6 \times 3} \mathbf N(x,y,z)\)称为几何矩阵
令\(\mathbf B_{6\times 12}(x,y) =\begin{bmatrix} \mathbf (B_1)_{6\times 3}& \mathbf (B_2)_{6\times 3} & \mathbf (B_3)_{6\times 3}& \mathbf (B_4)_{6\times 3}\end{bmatrix}\)
那么

\[\mathbf{B}_i=\begin{bmatrix}\partial\end{bmatrix} \begin{bmatrix}N_i & 0 & 0\\ 0 & N_i & 0\\ 0 & 0 & N_i\end{bmatrix}=\frac{1}{6V} \begin{bmatrix} b_i & 0 & 0\\ 0 & c_i & 0\\ 0 & 0 & d_i\\ c_i & b_i & 0\\0 & d_i & c_i\\ d_i & 0 & b_i\end{bmatrix},i=1,2,3,4 \]

• 应力场函数
  \(\mathbf {\sigma}(x,y,z)_{(6\times 1)}=\mathbf{ D_{(6\times 6)}} \cdot \epsilon=\mathbf{ D} \cdot (\mathbf B \cdot \mathbf q ^e)\)
  其中，\(\mathbf D\)称为弹性系数矩阵
  令\(\mathbf S=\mathbf D \cdot \mathbf B\)，并称\(\mathbf S\)为应力函数矩阵

4. 最小势能原理
   \(\prod ^e=U-W=\frac{1}{2} \int _{\Omega ^e} \sigma _{ij} \epsilon _{ij} d \Omega-[\int _{\Omega ^e} \bar b_i u _i d \Omega+\int _{S_p ^e} \bar p_i u _i dA]\)
   \(=\frac{1}{2} \mathbf q^{eT}(\int _{\Omega ^e} \mathbf B^{T}\mathbf D \mathbf B d \Omega)\mathbf q^{e}-(\int _{\Omega ^e} \mathbf N^{T} \mathbf{\bar b} d \Omega+\int _{S_p ^e} \mathbf N^{T} \mathbf{\bar p} dA)^T \mathbf q^{e}\)
   \(=\frac{1}{2} \mathbf q^{eT} \mathbf K ^e\mathbf q^{e}-\mathbf P^{eT} \mathbf q^{e}\)
   其中，单元刚度矩阵\(\mathbf K ^e_{(12\times 12)}=\int _{\Omega ^e} \mathbf B^{T}\mathbf D_{(6\times 6)} \mathbf B d \Omega\)

单元节点载荷 \(\mathbf P^{eT} =\int _{\Omega ^e} \mathbf N^{T} \mathbf{\bar b} d \Omega+\int _{S_p ^e} \mathbf N^{T} \mathbf{\bar p} dA\)

应用最小势能原理，\(\frac{\partial \prod(\mathbf q^{e}) }{\partial\mathbf q^{e}}=0\)
则获得刚度方程，\(\mathbf K ^e \cdot \mathbf q^{e}=\mathbf P^{e}\)

5. 空间4节点三角锥（四面体）单元性质
   • 单元内的应力和应变均为常数，又称为常应变长应力体元，与CST性质相同
     
     
   • 没有针对节点位移的坐标变换
   • 对于应变梯度较大的区域，单元划分应适当加密，否则导致较大误差

4节点四面体单元的matlab实现

三个函数

• 刚度矩阵Tetrahedron3D4Node_Stiffness(E,NU,xi,yi,zi,xj,yj,zj,xm,ym,zm,xn,yn,zn）：该函数计算单元的刚度矩阵，输入弹性模量E，泊松比NU，4个节点i、j、m、n的坐标xi,yi,zi,xj,yj,zj,xm,ym,zm,xn,yn,zn，输出单元刚度矩阵k(12x12)。
• 单元组装Tetrahedron3D4Node_Assembly(Kk,k,ij.m,n)：该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、j、m、n，输出整体度矩阵KK
• 单元应力Tetrahedron3D4Node_Stress(E,NU,xi,yi,zi,xj,yj,zj,xm,ym,zm,xn,yn,zn,u)：该函数计算单元的应力，输入弹性模量E，泊松比NU，4个节点i、j、m、n的坐标xi,yi,zi,xj,yj,zj,xm,ym,zm,xn,yn,zn，单元的位移列阵u(12x1），输出单元的应力stress，由于它为常应力单元，则单元的应力分量为Sx,Sy,Sz,Sxy,Syz,Szx

function k=Tetrahedron3D4Node_Stiffness(E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4)
%该函数计算单元的刚度矩阵%输入弹性模量E，泊松比NU
%输入4个节点i、j、m、n的坐标xi,yi,zi,xj.yj.zj.xm,ym,zm,xn,yn,zn
%输出单元刚度矩阵k(12X12）

%以下数组与书中的对应关系
%betai - bi
%gammai - ci
%deltai - di
%i=1:4

xyz = [1 x1 y1 z1 ; 1 x2 y2 z2 ; 1 x3 y3 z3 ; 1 x4 y4 z4]; V = det(xyz)/6;
mbeta1 = [1 y2 z2;1 y3 z3;1 y4 z4]; 
mbeta2 = [1 y1 z1;1 y3 z3;1 y4 z4]; 
mbeta3 = [1 y1 zl;1 y2 z2;1 y4 z4]; 
mbeta4 = [1 y1 z1;1 y2 z2;1 y3 z3];

mgamma1 = [1 x2 z2;1 x3 z3;1 x4 z4]; 
mgamma2 = [1 x1 z1;1 x3 z3;1 x4 z4]; 
mgamma3 = [1 x1 z1;1 x2 z2;1 x4 z4]; 
mgamma4 = [1 x1 z1;1 x2 z2;1 x3 z3]; 

mdelta1 = [1 x2 y2;1 x3 y3;1 x4 y4]; 
mdelta2 = [1 x1 y1;1 x3 y3;1 x4 y4]; 
mdelta3 = [1 x1 y1;1 x2 y2;1 x4 y4];
mdelta4 = [1 xl y1;1 x2 y2;1 x3 y3];

betal = -1*det(mbetal); 
beta2 = det(mbeta2); 
beta3 = -1*det(mbeta3);
beta4 = det(mbeta4);

gammal=det(mgamma1); 
gamma2=-1*det(mgamma2); 
gamma3 = det(mgamma3);
gamma4=-1*det(mgamma4); 

deltal = -1*det(mdeltal);
delta2 = det(mdelta2); 
delta3 = -1*det(mdelta3);
delta4 = det(mdelta4);

B1=[beta1 0 0;0 gamma1 0;0 0 deltal;
gammal betal 0;0 deltal gammal;deltal 0 betal];
B2=[beta2 0 0;0 gamma2 0 0;0 0 delta2;
gamma2 beta2 0;0 delta2 gamma2;delta2 0 beta2];
B3=[beta3 0 0;0 gamma3 0;0 0 delta3;
gamma3 beta3 0;0 delta3 gamma3;delta3 0 beta3];
B4=[beta4 00 ;0 gamma4 0;0 0 delta4;
gamma4 beta4 0;0 delta4 gamma4;delta4 0 beta4);

B = [B1 B2 B3 B4]/(6*V);
D = (E/((1+NU)*(1-2*NU)))*[1-NU NU NU 0 0 0;NU 1-NU NU 0 0 0;NU NU 1-NU 0 0 0;
0 0 0 (1-2*NU)/2 0 0;0 0 0 0 (1-2*NU)/2 0;0 0 0 0 0 (1-2*NU)/2]
k =abs(V)*B'*D*B;

function y=Tetrahedron3D4Node_Assembly(KK,k,i,j.m,n)
%该函数进行单元刚度矩阵的组装
%输入单元刚度矩阵k
%输入单元的节点编号i、j、m、n
%输出整体刚度矩阵KK

DOF=[3*i-2:3*i,3*j-2:3*j,3*m-2:3*m,3*n-2:3*n]
for n1=1:12 
	for n2=1:12
		KK(DOF(n1),DOF(n2))= KK(DOF(n1),DOF(n2))+k(n1,n2);
	end 
end 
y=KK;

function y=Tetrahedron3D4Node_Stress(E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,u)
%该函数计算单元的应力
%输入弹性模量E，泊松比NU
%输入4个节点i、j、m、n的坐标xi,yi,zi,xi.yj.zj,xm,ym,zm,xn,yn,zn
%输入单元的位移列阵u(12X1)
%输出单元的应力stress(6X1）
%由于它为常应力单元，应力分量为Sx,Sy,Sz,Sxy,Syz,Szx

xyz=[1 x1 y1 z1;1 x2 y2 z2;1 x3 y3 z3;1 x4 y4 z4];
V = det(xyz)/6;
mbetal = [1 y2 z2;1 y3 z3;1 y4 z4];
mbeta2 = [1 y1 z1;1 y3 z3;1 y4 z4]; 
mbeta3 = [1 y1 z1;1 y2 z2;1 y4 z4];
mbeta4 = [1 y1 z1;1 y2 z2;1 y3 z3];

mgamma1 = [1 x2 z2;1 x3 z3;1 x4 z4]; 
mgamma2 = [1 x1 z1;1 x3 z3;1 x4 z4]; 
mgamma3 = [1 x1 z1;1 x2 z2;1 x4 z4];
mgamma4 = [1 x1 z1;1 x2 z2;1 x3 z3]; 

mdelta1 = [1 x2 y2;1 x3 y3;1 x4 y4];
mdelta2 = [1 x1 y1;1 x3 y3;1 x4 y4]; 
mdelta3 = [1 x1 y1;1 x2 y2;1 x4 y4]; 
mdelta4 = [1 x1 y1;1 x2 y2;1 x3 y3];

betal = -1*det(mbetal);
beta2 = det(mbeta2); 
beta3 =-1*det(mbeta3); 
beta4 = det(mbeta4);

gammal = det(mgammal); 
gamma2 = -1*det(mgamma2); 
gamma3 = det(mgamma3); 
gamma4 = -1*det(mgamma4); 

deltal = -1*det(mdelta1);
delta2 = det(mdelta2); 
delta3 = -1*det(mdelta3); 
delta4 = det(mdelta4);

B1 = [beta1 0 0;0 gamma1 0;0 0 delta1;
gammal betal 0;0 deltal gammal;deltal 0 betal];
B2 = [beta2 0 0;0 gamma2 0;0 0 delta2;
gamma2 beta2 0;0 delta2 gamma2;delta2 0 beta2];
B3 = [beta3 0 0;0 gamma3 0;0 0 delta3;
gamma3 beta3 0;0 delta3 gamma3;delta3 0 beta3];
B4 = [beta4 0 0;0 gamma4 0;0 0 delta4;
gamma4 beta4 0;0 delta4 gamma4;delta4 0 beta4];
B = [B1 B2 B3 B4]/(6*V);
D = (E/(1+NU)*(1-2*NU)))*[1-NU NU NU 0 0 0;NU 1-NU NU 0 0 0;NU NU 1-NU 0 0 0;0 0 0 (1-2*NU)/2 0 0;0 0 0 0 (1-2*NU)/2 0;0 0 0 0 0 (1-2*NU)/2]
y=D*B*u;

空间8节点正六面体单元及MATLAB编程

每个节点具有3个自由度，共有24个DOF，整体坐标系建模

1. 节点描述

• 节点坐标：\((x_i,y_i,z_i),i=1,2,3,4,…,8\)
• 节点位移矩阵（8*\1）：
  \(\mathbf q^e=\begin{bmatrix} u_1 & v_1 & w_1 & u_2 & v_2 & w_2 & u_3 & v_3 & w_3 & … & u_8 & v_8 & w_8\end{bmatrix}^T\)
• 节点力矩阵（24*\1）：
  \(\mathbf P ^e=\begin{bmatrix} P_{x_1} & P_{y_1} & P_{z_1}& P_{x_{2}}& P_{y_2}& P_{z_2} & P_{x_3} & P_{y_{3}} & P_{z_3} & … & P_{x_8} & P_{y_{8}} & P_{z_8} \end{bmatrix}^T\)

2. 单元上的场描述

• 位移场函数：所有节点间使用相同函数插值，(低阶到高阶，唯一确定)
  \(u_i(x_i,y_i,z_i)=a_0+ a_1x_i+ a_2y_i+ a_3z_i+a_4x_iy_i+a_5y_iz_i+a_6z_ix_i+a_7x_iy_iz_i\)
  \(v_i(x_i,y_i,z_i)=b_0+ b_1x_i+ b_2y_i+ b_3z_i+b_4x_iy_i+b_5y_iz_i+b_6z_ix_i+b_7x_iy_iz_i\)
  \(w_i(x_i,y_i,z_i)=c_0+ c_1x_i+ c_2y_i+ c_3z_i+c_4x_iy_i+c_5y_iz_i+c_6z_ix_i+c_7x_iy_iz_i\)

  < 共有24个线性方程 >

  由\((u_i,v_i,w_i)\)可以确定系数\(( a_i, b_i, c_i)\)的值，其中\(i=1,2,3,4,…,8\)，得到系数，重新排列位移场
  \(u(x,y,z)=N_1(x,y,z)u_1+N_2u_2+N_3u_3+…+N_8u_8\) \(v(x,y,z)=N_1(x,y,z)v_1+N_3v_2+N_3v_3+…+N_8v_8\)
  \(w(x,y,z)=N_1(x,y,z)w_1+N_2w_2+N_3w_3+…+N_8w_8\)

\(N_i(x,y,z)\)的具体表达式可以由拉格朗日插值直接得出，自己去算

位移场矩阵形式：\(\mathbf u_{3\times 1}(x,y,z)=\mathbf N_{3\times 24}(x,y,z) \cdot \mathbf q ^e\)

其中\(\mathbf N(x,y,z)=\begin{bmatrix} N_1 & 0 & 0 & N_2 & 0 & 0 & … & N_8 & 0 & 0\\ 0 & N_1 & 0 & 0 & N_2 & 0 & … & 0 & N_8 & 0 \\ 0 & 0 & N_1 & 0 & 0 & N_2 & … & 0 & 0 & N_8 \end{bmatrix}\)

• 应变场函数

  \(\mathbf{\epsilon}(x,y,z)=\begin{bmatrix}\partial \end{bmatrix}_{(6 \times 3)} \mathbf{u}(x,y)=\begin{bmatrix}\partial \end{bmatrix}_{(6 \times 3)} \mathbf N(x,y) \cdot \mathbf q ^e=\mathbf B_{(6\times 24)}(x,y) \cdot \mathbf q ^e\)
  其中，\(\begin{bmatrix}\partial\end{bmatrix}=\begin{bmatrix} \frac{\partial}{\partial x} & 0 & 0\\ 0 & \frac{\partial}{\partial y} & 0\\0 & 0 & \frac{\partial}{\partial z}\\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0\\ 0 & \frac{\partial}{\partial z}& \frac{\partial}{\partial y}\\ \frac{\partial}{\partial z} & 0 & \frac{\partial}{\partial x}\end{bmatrix}\)，\(\mathbf B_{6\times 24}(x,y,z)=\begin{bmatrix}\partial \end{bmatrix}_{6 \times 3} \mathbf N(x,y,z)\)称为几何矩阵
  令\(\mathbf B_{6\times 24}(x,y) =\begin{bmatrix} \mathbf (B_1)_{6\times 3}& \mathbf (B_2)_{6\times 3} & … & \mathbf (B_8)_{6\times 3}\end{bmatrix}\)

  那么

  \[\mathbf{B}_i=\begin{bmatrix}\partial\end{bmatrix} \begin{bmatrix}N_i & 0 & 0\\ 0 & N_i & 0\\ 0 & 0 & N_i\end{bmatrix},i=1,2,3,4,…,8 \]

• 应力场函数
  \(\mathbf {\sigma}(x,y,z)_{(6\times 1)}=\mathbf{ D_{(6\times 6)}} \cdot \epsilon=\mathbf{ D} \cdot (\mathbf B \cdot \mathbf q ^e)=\mathbf S{_{(6\times 24)}} \cdot \mathbf q ^e\)
  其中，\(\mathbf D\)称为弹性系数矩阵
  令\(\mathbf S=\mathbf D \cdot \mathbf B\)，并称\(\mathbf S\)为应力函数矩阵

4. 最小势能原理
   \(\prod ^e=U-W=\frac{1}{2} \int _{\Omega ^e} \sigma _{ij} \epsilon _{ij} d \Omega-[\int _{\Omega ^e} \bar b_i u _i d \Omega+\int _{S_p ^e} \bar p_i u _i dA]\)
   \(=\frac{1}{2} \mathbf q^{eT}(\int _{\Omega ^e} \mathbf B^{T}\mathbf D \mathbf B d \Omega)\mathbf q^{e}-(\int _{\Omega ^e} \mathbf N^{T} \mathbf{\bar b} d \Omega+\int _{S_p ^e} \mathbf N^{T} \mathbf{\bar p} dA)^T \mathbf q^{e}\)
   \(=\frac{1}{2} \mathbf q^{eT} \mathbf K ^e\mathbf q^{e}-\mathbf P^{eT} \mathbf q^{e}\)

其中，单元刚度矩阵\(\mathbf K ^e_{(24\times 24)}=\int _{\Omega ^e} \mathbf B^{T}\mathbf D_{(6\times 6)} \mathbf B d \Omega\)

单元节点载荷 \(\mathbf P^{eT} =\int _{\Omega ^e} \mathbf N^{T} \mathbf{\bar b} d \Omega+\int _{S_p ^e} \mathbf N^{T} \mathbf{\bar p} dA\)

应用最小势能原理，\(\frac{\partial \prod(\mathbf q^{e}) }{\partial\mathbf q^{e}}=0\)
则获得刚度方程，\(\mathbf K ^e \cdot \mathbf q^{e}=\mathbf P^{e}\)

5. 空间4节点三角锥（四面体）单元性质

• 位移场为完全线性+交叉项函数
  

• 应变场不完全线性函数
  

• 应力场为不完全线性函数：\(\mathbf {\sigma}=\mathbf{ D_{(6\times 6)}} \cdot \epsilon\)

• (1) 位移场为完全线性+交叉项函数，在x，y，z方向呈线性变化

• (2) 应变场及应力场为不完全线性函数

• (3) 空间8节点正六面体单元的比4节点四面体单元的精度高；与平面4节点矩形单元的性质相同

• (4) 正六面体单元的几何形状还应变换为任意六面体的形状，以增强几何的适应性

8节点正六面体单元的matlab实现

三个函数

• Hexahedral3D8Node_Stiffness(E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8)
  • 该函数计算单元的刚度矩阵，输入弹性模量E，泊松比NU，8个节点的坐标x1、y1、z1、x2、y2、z2、x3、y3、z3、x4、y4、z4、x5、y5、z5、x6、y6、z6、x7、y7、z7、x8、y8、z8，输出单元刚度矩阵k(24X24)
• Hexahedral3D8Node_Assembly(kk,k,i,j,l,m,n,o,p,q)
  • 该函数进行单元刚度矩阵的组装，输入单元刚度矩阵k，单元的节点编号i、j、I、m、n、o、p、q，输出整体刚度矩阵kk。
• Hexahedral3D8Node_Stress(E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,u)
  • 该函数计算单元中心点的应力，输入弹性模量E，泊松比NU，8个节点的坐标x1、y1、z1、x2、y2、z2、x3、y3、z3、x4、y4、z4、x5、y5、z5、x6、y6、z6、x7、y7、z7、x8、y8、z8，单元的位移列阵u(6X1)，输出单元中心的应力stress(6x1），则单元的应力分量为Sx,Sy,Sz,Sxy,Syz,Szx

刚度矩阵函数

矩阵太长的有换行，注意使用时删除

function k=Hexahedral3D8Node_Stiffness(E,NU,x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,
x6,y6,z6,x7,y7,z7,x8,y8,z8)
%该函数计算单元的刚度矩阵
%输入弹性模量E，泊松比NU
%输入8个节点的坐标x1、y1、z1、x2、y2、z2、x3、y3、z3、x4、y4、z4、
%输入x5、y5、z5、x6、y6、z6、x7、y7、z7、x8、y8、z8
%输出单元刚度矩阵k(24X24)

syms s t n;   %定义形状函数N

%定义局部坐标系
N1=(1+s)*(1-t)*(1-n)/8; 
N2=(1+s)*(1+t)*(1-n)/8;
N3=(1-s)*(1+t)*(1-n)/8; 
N4=(1-s)*(1-t)*(1-n)/8; 
N5=(1+s)*(1-t)*(1+n)/8;
N6=(1+s)*(1+t)*(1+n)/8;
N7=(1-s)*(1+t)*(1+n)/8;
N8=(1-s)*(1-t)*(1+n)/8;

%定义坐标变换
X=N1*×1+N2*×2+N3*×3+N4*×4+N5*×5+N6*×6+N7*×7+N8*×8; y=N1*y1+N2*y2+N3*y3+N4*y4+N5*y5+N6*y6+N7*y7+N8*y8;
z=N1*z1+N2*z2+N3*z3+N4*z4+N5*z5+N6*z6+N7*z7+N8*z8;

%定义雅可比矩阵
J=[diff(x,s),diff(y,s),diff(z,s);
diff(x,t),diff(y,t),diff(z,t);
diff(x,n),diff(y.n),diff(z,n)];

Jdet=det(J); 
J %打印J
Jdet %打印Jdet

%定义B矩阵的系数
a=diff(y,t)*diff(z,n)-diff(z,t)*diff(y,n);
b=diff(y,s)*diff(z,n)-diff(z,s)*diff(y,n);
c=diff(y,s)*diff(z,t)-diff(z,s)*diff(y,t); 
d=diff(x,t)*diff(z,n)-diff(z,t)*diff(x,n);
e=diff(x,s)*diff(z.n)-diff(z,s)*diff(x,n); 
f=diff(x,s)*diff(z,t)-diff(z,s)*diff(x,t); 
g=diff(x,t)*diff(y.n)-diff(y,t)*diff(x,n);
h=diff(x,s)*diff(y.n)-diff(y,s)*diff(x,n);
I=diff(x,s)*diff(y.t)-diff(y,s)*diff(x,t);

%通过循环计算各个矩阵
Ns=[N1,N2,N3,N4,N5,N6,N7,N8];
Bs=sym(zeros(6,3,8));
for i=1:8
	Bs(:,:,i)=[a*diff(Ns(i),s)-b*diff(Ns(i),t)+c*diff(Ns(i),n),0,0;
	0,-d*diff(Ns(i),s)+e*diff(Ns(i),t)-f*diff(Ns(0),n),0;
	0,0,g*diff(Ns(0),s)-h*diff(Ns(0),t)+I*diff(Ns(i),n); 
	-d*diff(Ns(i),s)+e*diff(Ns(i),t)-f*diff(Ns(i),n),a*diff(Ns(i),s)-b*diff(Ns(i),t)+c*diff(Ns(i),n),0;
O,g*diff(Ns(0),s)-h*diff(Ns(i),t)+I*diff(Ns(i),n),-d*diff(Ns(i),s)+e*diff(Ns(i),t)-f*diff(Ns(i),n);
g*diff(Ns(i),s)-h*diff(Ns(i),t)+I*diff(Ns(i),n),O,a*diff(Ns(i),s)-b*diff(Ns(i),t)+c*diff(Ns(i),n)]/Jdet;
end

%计算B矩阵
B=[Bs(:,:,1),Bs(:,:,2),Bs(:,:,3),Bs(:,:,4),Bs(:,:,5),Bs(:,:,6),Bs(:,:,7),Bs(:,:,8)]
B

%计算D矩阵
D=(E/((1+NU)*(1-2*NU)))*[1-NU,NU,NU,0,0,0;NU,1-NU,NU,0,0,0;NU,NU,1-NU,0,0,0;
0,0,0,0.5-NU,0,0;0,0,0,0,0.5-NU,0;0,0,0,0,0,0.5-NU); 
D

%计算单元刚度矩阵k
BD=Jdet*transpose(B)*D*B;
z=(int(int(int(BD,n,-1,1),t,-1,1),s,-1,1));
z
k=double(z);

刚度矩阵组装函数

functionz=Hexahedral3D8Node_Assembly(KK,k,i,j,l,m,n,o,p,q)
%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k
%输入单元的节点编号i、j、l、m、n、o、p、q
%输出整体刚度矩阵KK

DOF=[3*i-2:3*i,3*j-2:3*j,3*l-2:3*|,3*m-2:3*m,3*n-2:3*n,3*o-2:3*o,3*p-2:3*p,3*q-2:3*q];
for n1=1:24 
	for n2=1:24
		KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);
	end 
end 
z=KK;

应力计算函数

functionstress=Hexahedral3D8NodeStress(E,NU,x1,y1,z1.x2.y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,u)
%该函数计算单元中心点的应力
%输入弹性模量E，泊松比NU
%输入8个节点的坐标x1、y1、z1、x2、y2、z2、x3、y3、z3、x4、y4、z4、
%输入x5、y5、z5、x6、y6、z6、x7、y7、z7、x8、y8、z8
%输入单元的位移列阵u（6X1）
%输出单元中心的应力stress(6X1),应力分量为Sx,Sy,Sz,Sxy.Syz,Szx

syms s t n;  %定义局部坐标系

%定义形状函数N
N1=(1+s)*(1-t)*(1-n)/8; 
N2=(1+s)*(1+t)*(1-n)/8;
N3=(1-s)*(1+t)*(1-n)/8; 
N4=(1-s)*(1-t)*(1-n)/8; 
N5=(1+s)*(1-t)*(1+n)/8; 
N6=(1+s)*(1+t)*(1+n)/8; 
N7=(1-s)*(1+t)*(1+n)/8; 
N8=(1-s)*(1-t)*(1+n)/8;

%定义坐标变换
x=N1*x1+N2*x2+N3*x3+N4*x4+N5*×5+N6*x6+N7*x7+N8*x8;
y=N1*y1+N2*y2+N3*y3+N4*y4+N5*y5+N6*y6+N7*y7+N8*y8;
z=N1*z1+N2*z2+N3*z3+N4*z4+N5*z5+N6*z6+N7*z7+N8*z8;

%定义雅可比矩阵
J=[diff(x,s),diff(y,s),diff(z,s);diff(x,t),diff(y.t),diff(z,t);diff(x,n),diff(y,n),diff(z,n)];
Jdet=det(J);

%定义B矩阵的系数
a=diff(y,t)*diff(z,n)-diff(z,t)*diff(y,n);
b=diff(y,s)*diff(z,n)-diff(z,s)*diff(y,n);
c=diff(y,s)*diff(z,t)-diff(z,s)*diff(y,t); 
d=diff(x,t)*diff(z,n)-diff(z,t)*diff(x,n);
e=diff(x,s)*diff(z.n)-diff(z,s)*diff(x,n); 
f=diff(x,s)*diff(z,t)-diff(z,s)*diff(x,t); 
g=diff(x,t)*diff(y.n)-diff(y,t)*diff(x,n);
h=diff(x,s)*diff(y.n)-diff(y,s)*diff(x,n);
I=diff(x,s)*diff(y.t)-diff(y,s)*diff(x,t);

%通过循环计算各个矩阵
Ns=[N1,N2,N3,N4,N5,N6,N7,N8];
Bs=sym(zeros(6,3,8));
for i=1:8
	Bs(:,:,i)=[a*diff(Ns(i),s)-b*diff(Ns(i),t)+c*diff(Ns(i),n),0,0;
	0,-d*diff(Ns(i),s)+e*diff(Ns(i),t)-f*diff(Ns(0),n),0;
	0,0,g*diff(Ns(0),s)-h*diff(Ns(0),t)+I*diff(Ns(i),n); 
	-d*diff(Ns(i),s)+e*diff(Ns(i),t)-f*diff(Ns(i),n),a*diff(Ns(i),s)-b*diff(Ns(i),t)+c*diff(Ns(i),n),0;
O,g*diff(Ns(0),s)-h*diff(Ns(i),t)+I*diff(Ns(i),n),-d*diff(Ns(i),s)+e*diff(Ns(i),t)-f*diff(Ns(i),n);
g*diff(Ns(i),s)-h*diff(Ns(i),t)+I*diff(Ns(i),n),O,a*diff(Ns(i),s)-b*diff(Ns(i),t)+c*diff(Ns(i),n)]/Jdet;
end

%计算B矩阵
B=[Bs(:,:,1),Bs(:,:,2),Bs(:,:,3),Bs(:,:,4),Bs(:,:,5),Bs(:,:,6),Bs(:,:,7),Bs(:,:,8)]

%计算D矩阵
D=(E/((1+NU)*(1-2*NU)))*[1-NU,NU,NU,0,0,0;NU,1-NU,NU,0,0,0;NU,NU,1-NU,0,0,0;
0,0,0,0.5-NU,0,0;0,0,0,0,0.5-NU,0;0,0,0,0,0,0.5-NU); 

%计算应力向量 
W=D*B*u;
%在单元的中心处计算应力 
wcent=subs(w,{s,t,n},{0,0,0}) 
stress=double(wcent)

参数单元的原理

具有曲边的单元划分：以平面问题的任意四边形单元为例，内部可以采用矩形单元，逼近曲边会产生锯齿

采用任意四边形单元可以较好地逼近曲边

• 方法一：直接构建
• 方法二：基于矩形单元来进行几何形状的映射（参数单元）

模型构建采用基准单元，坐标系\((\xi ,\eta)\)

实际计算时采用物理单元，坐标系\((x ,y)\)

通过几何形状映射（3个方面）来进行计算：参数单元

1. 几何坐标函数映射：\((x,y)\Longrightarrow (\xi ,\eta )\)

坐标映射

映射关系：采用多项式表达

\(\left\{\begin{matrix} x=a_0+a_1\xi+a_2\eta +a_3 \xi \eta\\y=b_0+b_1\xi+b_2\eta +b_3 \xi \eta\end{matrix}\right.\)

位移场

\(\left\{\begin{matrix} x(\xi, \eta)=\widetilde{N}_1(\xi, \eta)x_1+\widetilde{N}_2(\xi, \eta)x_2+\widetilde{N}_3(\xi, \eta)x_3+\widetilde{N}_4(\xi, \eta)x_4 \\y(\xi, \eta)=\widetilde{N}_1(\xi, \eta)y_1+\widetilde{N}_2(\xi, \eta)y_2+\widetilde{N}_3(\xi, \eta)y_3+\widetilde{N}_4(\xi, \eta)y_4\end{matrix}\right.\)

其中，\(\widetilde{N}_i=\frac{1}{4}(1+\xi _i \xi)(1+\eta _i \eta ),i=1,2,3,4\)

\((x,y)\)坐标系下的节点坐标值排列：

\(\mathbf {\widetilde{q}}=\begin{bmatrix} x_1 & y_1 & x_2 & y_2 & x_3 & y_3 & x_4 & y_4\end{bmatrix}^T\)

则坐标映射可写成：\((x,y)\Longrightarrow (\xi ,\eta )\)

\(\begin{bmatrix}x \\ y\end{bmatrix}=\begin{bmatrix}x(\xi,\eta) \\ y(\xi,\eta)\end{bmatrix}=\begin{bmatrix} \widetilde{N}_1 & 0 & \widetilde{N}_2 & 0 & \widetilde{N}_3 & 0 & \widetilde{N}_4 & 0\\ 0 & \widetilde{N}_1 & 0 & \widetilde{N}_2 & 0 & \widetilde{N}_3 & 0 & \widetilde{N}_4\end{bmatrix}=\mathbf {\widetilde{N}}_{(2\times 8)}\cdot \mathbf {\widetilde{q}}\)
2. 坐标的偏导数变换：\((\frac{\partial }{\partial x} ,\frac{\partial }{\partial y})\Longrightarrow (\frac{\partial }{\partial \xi} ,\frac{\partial }{\partial \eta })\)

\(\left\{\begin{matrix} \frac{\partial \Phi }{\partial \xi} =\frac{\partial \Phi }{\partial x}\frac{\partial x }{\partial \xi}+\frac{\partial \Phi }{\partial y}\frac{\partial y }{\partial \xi}\\ \frac{\partial \Phi }{\partial \eta} =\frac{\partial \Phi }{\partial x}\frac{\partial x }{\partial \eta}+\frac{\partial \Phi }{\partial y}\frac{\partial y }{\partial \eta}\end{matrix}\right.\)

省略目标函数\(\Phi\)，并交换位置

\(\left\{\begin{matrix} \frac{\partial }{\partial \xi} =\frac{\partial x }{\partial \xi}\frac{\partial }{\partial x}+\frac{\partial y }{\partial \xi}\frac{\partial }{\partial y}\\ \frac{\partial }{\partial \eta} =\frac{\partial x }{\partial \eta}\frac{\partial }{\partial x}+\frac{\partial y }{\partial \eta}\frac{\partial }{\partial y}\end{matrix}\right.\)

改写成矩阵形式

\(\begin{bmatrix} \frac{\partial }{\partial \xi} \\ \frac{\partial }{\partial \eta}\end{bmatrix} = \mathbf{J} \begin{bmatrix} \frac{\partial }{\partial x} \\ \frac{\partial }{\partial y}\end{bmatrix}\)

其中，Jacobi矩阵\(\mathbf{J}=\begin{bmatrix}\frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi}\\ \frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta}\end{bmatrix}\)，\(|\mathbf {J}|=\frac{\partial x}{\partial \xi} \frac{\partial y}{\partial \eta}- \frac{\partial y}{\partial \xi} \frac{\partial x}{\partial \eta}\)

求逆，即可完成偏导数的映射

\(\begin{bmatrix} \frac{\partial }{\partial x }\\ \frac{\partial }{\partial y }\end{bmatrix}={\mathbf {J}}^{-1}\begin{bmatrix} \frac{\partial }{\partial \xi }\\ \frac{\partial }{\partial \eta }\end{bmatrix}=\frac{1}{|\mathbf {J}|} \begin{bmatrix}\frac{\partial y}{\partial \eta} & -\frac{\partial y}{\partial \xi}\\ -\frac{\partial x}{\partial \eta} & \frac{\partial x}{\partial \xi}\end{bmatrix} \begin{bmatrix} \frac{\partial }{\partial \xi }\\ \frac{\partial }{\partial \eta }\end{bmatrix}\)

方程形式

\(\left\{\begin{matrix}\frac{\partial }{\partial x} =\frac{1}{|\mathbf {J}|}(\frac{\partial y}{\partial \eta}\frac{\partial }{\partial \xi}- \frac{\partial y}{\partial \xi}\frac{\partial }{\partial \eta}) \\\frac{\partial }{\partial y} =\frac{1}{|\mathbf {J}|}(\frac{\partial x}{\partial \eta}\frac{\partial }{\partial \xi}+ \frac{\partial x}{\partial \xi}\frac{\partial }{\partial \eta})\end{matrix}\right.\)

3. 面积和体积的变换：\(\int _{A^e}dxdy\Longrightarrow \int _{A^e}d\xi d\eta\)

物理坐标系\((x,y)\)中的单位矢量为\((\mathbf{i},\mathbf{j})\)

面积微元=坐标微元的叉乘，即\(dA=|d\xi \times d\eta|\)

\(\left\{\begin{matrix}d\mathbf{\xi} = \frac{\partial x}{\partial \xi} \partial \xi \cdot \mathbf{i}+ \frac{\partial y}{\partial \xi} \partial \xi \cdot \mathbf{j} \\d\mathbf{\eta} = \frac{\partial x}{\partial \eta} \partial \eta \cdot \mathbf{i}+ \frac{\partial y}{\partial \eta} \partial \eta \cdot \mathbf{j}\end{matrix}\right.\)

矩阵形式

\(dA=|\mathbf {J}|d\xi d\eta =\begin{bmatrix}\frac{\partial x}{\partial \xi} d\xi & \frac{\partial y}{\partial \xi}d\xi\\ \frac{\partial x}{\partial \eta}d\eta & \frac{\partial y}{\partial \eta}d\eta\end{bmatrix}\)

类似，三维问题中的体积微元=两坐标微元叉积点乘第三坐标微元，即

\(d\Omega=d\xi \cdot (d\xi \times d\zeta)\)

\(d\Omega =\begin{vmatrix} \frac{\partial x}{\partial \xi} d\xi & \frac{\partial y}{\partial \xi} d\xi & \frac{\partial z}{\partial \xi} d\xi\\ \frac{\partial x}{\partial \eta} d\eta & \frac{\partial y}{\partial \eta} d\eta & \frac{\partial z}{\partial \eta} d\eta\\\frac{\partial x}{\partial \zeta} d\zeta & \frac{\partial y}{\partial \zeta} d\zeta & \frac{\partial z}{\partial \zeta} d\zeta\end{vmatrix}=|\mathbf{J} |d\xi d\eta d\zeta\)

4. 参数单元（平面单元）

三个方面的坐标变换

• \((x,y)\Longrightarrow (\xi ,\eta )\)
• \((\frac{\partial }{\partial x} ,\frac{\partial }{\partial y})\Longrightarrow (\frac{\partial }{\partial \xi} ,\frac{\partial }{\partial \eta })\)
• \(\int _{A^e}dxdy\Longrightarrow \int _{A^e}d\xi d\eta\)

\((x,y)\)坐标系下的单元几何矩阵的变换

\(\mathbf{B}(x,y,\frac{\partial }{\partial x} ,\frac{\partial }{\partial y})=\begin{bmatrix} \frac{\partial }{\partial x} & 0\\ 0 & \frac{\partial }{\partial y}\\ \frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{bmatrix}\mathbf{B}(x,y)=\mathbf{B^*}(\xi,\eta ,\frac{\partial }{\partial \xi} ,\frac{\partial }{\partial \eta})\)

\((x,y)\)坐标系下的单元刚度矩阵的变换

\(\mathbf{K}^e_{(xy)} =\int_{A^e} \mathbf{B}^T(x,y,\frac{\partial }{\partial x} ,\frac{\partial }{\partial y})\cdot \mathbf{D} \cdot \mathbf{B}(x,y,\frac{\partial }{\partial x} ,\frac{\partial }{\partial y})dA\cdot t\)
\(=\int_{A^e} \mathbf{B}^{*T}(x,y,\frac{\partial }{\partial x} ,\frac{\partial }{\partial y})\cdot \mathbf{D} \cdot \mathbf{B}^{*}(x,y,\frac{\partial }{\partial x} ,\frac{\partial }{\partial y}) |\mathbf {J}|d\xi d\eta \cdot t\)

将\((x,y)\)坐标系下全部转换到\((\xi ,\eta )\)坐标系下进行计算

5. 参数单元的三种类型

几何形状插值函数的阶次与位移插值函数的阶次对比

数值积分

数值积分举例——Gauss积分

\(I=\int _{-1}^1 f(\xi) d\xi \approx \sum _{k=1}^n A_kf(\xi _k)\)

• 积分权系数\(A_k\)
• 积分点位置\(\xi _k\)

最佳数值积分基本思想，构造多项式\(\varphi(\xi _i)\)去等值\(f(\xi _i)\)，即

\(\varphi(\xi _i)=f(\xi _i);(i=1,2,3,…,n)\)

\(i\)取值足够大时，则有近似代替

\(I=\int _{-1}^1 f(\xi) d\xi \approx \int _{-1}^1 \varphi (\xi) d\xi\)

技术关键：如何构造最好的多项式\(\varphi(\xi _i)\)去逼近原函数，这里采用Gauss数值积分

Gauss数值积分

• 1点Gauss数值积分公式（\(n=1\)）

  \(\int _{-1}^1 f(\xi) d\xi \approx \sum _{k=1}^1 A_kf(\xi _k)=A_1f(\xi _1)\)

  使用两种多项式（常数，一次）代替\(f(\xi )\)

  • 令\(f(\xi )=a_0\)时，有 \(\int _{-1}^1 a_0 d\xi =A_1\cdot a_0\)
  • 令\(f(\xi )=a_1\xi\)时，有 \(\int _{-1}^1 a_1\xi d\xi =A_1\cdot( a_1\xi _1)\)
    解得，\(A_1=2，\xi _1=0\)

  最终的1点Gauss积分（体形积分）为：\(\int _{-1}^1 f(\xi) d\xi \approx =2f(0)\)
  

• 2点Guass积分公式（\(n=2\)）

  \(\int _{-1}^1 f(\xi) d\xi \approx A_1f(\xi _1)+A_2f(\xi _2)\)

  4个系数确定有两种方法

  • 构造正交多项式

  • 直接法（以下采用这方法）

    选择四个多项式(\(1,\xi ,\xi ^2,\xi ^3\))去代替\(f(\xi )\)，得到四个方程

    \(\left.\begin{matrix}2=A_1+A_2 \\0=A_1\xi _1+A_2\xi _2 \\\frac{2}{3}=A_1\xi _1^2+A_2\xi _2^2 \\0=A_1\xi _1^3+A_2\xi _2^3\end{matrix}\right\}\)

    解得

    \(\left.\begin{matrix}\xi _1=-\frac{1}{\sqrt{3} } \\\xi _2=\frac{1}{\sqrt{3} } \\A_1=1 \\A_2=1\end{matrix}\right\}\)

    所以，\(\int _{-1}^1 f(\xi) d\xi \approx f(-\frac{1}{\sqrt{3} })+f(\frac{1}{\sqrt{3} })\)

• n点Gauss积分

  \(\int _{-1}^1 f(\xi) d\xi \approx A_1f(\xi _1)+A_2f(\xi _2)+…+A_nf(\xi _n)\)

  采用直接设定多项式的方法求解过程较为复杂，采用Legendre多项式求取系数，查询Gauss数值积分手册

  

刚度矩阵的数值积分

在\((\xi ,\eta )\)坐标系下进行计算刚度矩阵

\(\mathbf{K}^e_{(xy)} =\int_{A^e} \mathbf{B}^T(x,y,\frac{\partial }{\partial x} ,\frac{\partial }{\partial y})\cdot \mathbf{D} \cdot \mathbf{B}(x,y,\frac{\partial }{\partial x} ,\frac{\partial }{\partial y})dA\cdot t\)
\(=\int_{A^e} \mathbf{B}^{*T}(x,y,\frac{\partial }{\partial x} ,\frac{\partial }{\partial y})\cdot \mathbf{D} \cdot \mathbf{B}^{*}(x,y,\frac{\partial }{\partial x} ,\frac{\partial }{\partial y}) |\mathbf {J}|d\xi d\eta \cdot t\)

2D平面4节点等参元刚度矩阵元素

\(k^e_{(xy)ij}=\int _{-1}^1 \int _{-1}^1 \frac{1}{A_0+B_0\xi +C_0 \eta}[(A_{\alpha i}+B_{\alpha i}\xi +C_{\alpha i} \eta)(A_{\beta j}+B_{\beta j}\xi +C_{\beta j} \eta)]d\xi d\eta \cdot t\)
对其需要进行数值积分，这是二重积分问题，采用Gauss数值积分

\(I=\int _{-1}^1 \int _{-1}^1 f(\xi ,\eta)d\xi d\eta=\sum _{j=1}^n \sum_{i=1}^n A_jA_i f(\xi _j,\eta _i)=\sum _{i,j=1}^n A_{ij} f(\xi _j,\eta _i)\)

其中，\(A_{ij}=A_iA_j\)

对于3D问题类似，

\(I=\int _{-1}^1 \int _{-1}^1 \int _{-1}^1 f(\xi ,\eta ,\zeta)d\xi d\eta d\zeta=\sum _{m=1}^n \sum _{j=1}^n \sum_{i=1}^n A_jA_i f(\xi _j,\eta _i,\zeta)=\sum _{i,j,m=1}^n A_{ijm} f(\xi _j,\eta _i,\zeta)\)

其中，\(A_{ijm}=A_iA_jA_m\)

典型空间问题的MATLAB编程

如图所示的一个空间块体，左端固定在面上，右端部受两个集中力F作用。试用MATLAB程序计算计算各个节点位移、支座反力以及单元的应力。

结构的参数

• \(F =100000N\)
• \(E=1×10 ^{10}Pa\)
• \(\mu=0.25\)
• \(t =0.2m\)

有限元模型

问题求解

1. 结构的离散化与编号

单元号	节点号
1	1 4 2 6
2	1 4 3 7
3	6 7 5 1
4	6 7 8 4
5	1 4 6 7
总节点位移矩阵（24*\1）：
\(\mathbf q=\begin{bmatrix} u_1 & v_1 & w_1 & u_2 & v_2 & w_2 & u_3 & v_3 & w_3 & … & u_8 & v_8 & w_8\end{bmatrix}^T\)

8个节点外载荷：
\(\mathbf F=\begin{bmatrix} 0 & 0 & \mathbf F_3^T & \mathbf F_4^T & 0 & 0 & \mathbf F_7^T & \mathbf F_8^T \end{bmatrix}^T\)

其中，\(\mathbf F_3^T =\mathbf F_4^T =\begin{bmatrix} 0 & 0 & 0 \end{bmatrix}\)，\(\mathbf F_7^T =\mathbf F_8^T =\begin{bmatrix} 0 & 0 & -10^5N \end{bmatrix}\)

支反力，墙面对节点1、3、5、6存在支反力：
\(\mathbf R=\begin{bmatrix} \mathbf R_1^T & \mathbf R_2^T &0 & 0 & \mathbf R_5^T & \mathbf R_6^T & 0 & 0 \end{bmatrix}^T\)

其中，\(\mathbf R_1^T =\begin{bmatrix} R_{1x} & R_{1y} & R_{1z} \end{bmatrix}\)，\(\mathbf R_2^T =\begin{bmatrix} R_{2x} & R_{2y} & R_{2z} \end{bmatrix}\)，

\(\mathbf R_5^T =\begin{bmatrix} R_{5x} & R_{5y} & R_{5z} \end{bmatrix}\)，\(\mathbf R_6^T =\begin{bmatrix} R_{6x} & R_{6y} & R_{6z} \end{bmatrix}\)

总节点力\(\mathbf P=\mathbf F+\mathbf R=\begin{bmatrix} \mathbf R_1^T & \mathbf R_2^T & \mathbf F_3^T & \mathbf F_4^T & \mathbf R_5^T & \mathbf R_6^T & \mathbf F_7^T & \mathbf F_8^T \end{bmatrix}^T\)

matlab编程实现

2. 计算各单元的刚度矩阵

首先在MATLAB环境下，输入弹性模量E、泊松比NU，然后针对单元1到单元5，分别调用5次函数 Tetrahedron3D4Node_Stiffness，就可以得到单元的刚度矩阵k1(6x6）～k5(6×6)。

MATLAB窗口的输入、输出情况

>> E=1e10; >> NU=0.25;
>> k1=Tetrahedron3D4Node_Stiffness(E,NU,0,0,0,0.2,0.8,0,0.2,0,0,0.2,0,0.6) 
>> k2=Tetrahedron3D4Node_Stiffness(E,NU,0,0,0,0.2,0.8,0,0,0.8,0,0,0.8,0.6); 
>>  k3=Tetrahedron3D4Node_Stiffness(E,NU,0.2,0,0.6,0,0.8,0.6,0,0,0.6,0,0,0);
>>  k4=Tetrahedron3D4Node_Stiffness(E,NU,0.2,0,0.6,0,0.8,0.6,0.2,0.8,0.6,0.2,0.8,0);
>>  k5=Tetrahedron3D4Node_Stiffness(E,NU,0,0,0,0.2,0.8,0,0.2,0,0.6,0,0.8,0.6);

3. 建立整体刚度矩阵

KK清零，然后5次调用函数Tetrahedron3D4Node_Assembly进行刚度矩阵的组装。

MATLAB窗口的输入、输出情况

>> Kk=zeros(24)
>> KK = Tetrahedron3D4Node_Assembly(KK,k1,1,4,2,6); 
>> KK = Tetrahedron3D4Node_Assembly(KK,k2,1,4,3,7); 
>> KK = Tetrahedron3D4Node_Assembly(KK,k3,6,7,5,1); 
>> KK = Tetrahedron3D4Node_Assembly(KK,k4,6,7,8,4); 
>> KK = Tetrahedron3D4Node_Assembly(KK,k5,1,4,6,7);

4. 边界条件的处理及刚度方程求解

节点1、2、5、6上三个方向的位移将为零，因此，将针对节点3、4、7、8的位移进行求解

节点1、2、5、6的位移将对应KK矩阵中的第1至6行，第13至18行和第1至6列，第13至18列，需从KK(24×24)中提出，置给k，然后生成对应的载荷列p，用高斯消去法进行求解

注意：MATLAB中的反斜线符号“”就是采用高斯消去法

>> k=KK([7:12,19:24],[7:12,19:24]); 
>> p=[0,0,0,0,0,0,0,0,-1e5,0,0,-1e5]';
>> u=k\p
u=1.0e-003 *
0.1249 -0.0485 -0.4024 0.1343 -0.0715 -0.4031
0.1314 0.0858 -0.4460 0.1353 0.0681 -0.4742

空间块体的节点位移计算结果(m)

5. 支反力的计算

在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力：先将上面得到的位移结果与位移边界条件的节点位移进行组合（注意位置关系），可以得到整体的位移列U(24×1)再代回原整体刚度方程，计算出所有的节点力P(24×1），接按对应关系就可以找到对应的支反力

>> U=[zeros(6,1);u([1:6]);zeros(6,1);u(7:12)]; 
>> P=KK*U
P = 1.0e+005 *
0.3372 1.3774 0.1904 -0.4202 1.2892 0.4984
-0.0000 0.0000 0.0000 -0.0000 -0.0000-0.0000
-0.4745 -1.3774 0.5604 0.5575 -1.2892 0.7509
-0.0000 -0.0000 -1.0000 -0.0000 0.0000 -1.0000

空间块体的支反力计算结果

6. 各单元的应力计算

先从整体位移列阵U(24×1)中提取出单元的位移列阵，然后，调用计算单元应力的函数 Tetrahedron3D4Node_Stress，就可以得到各个单元的应力分量。

MATLAB窗口的输入、输出情况

空间块体的各个单元应力分量的计算结果

典型空间问题的ANSYS分析实例

使用ansys的APDL语言处理上述问题

APDL可以实现-模型的参数化

• 获取ANSYS数据库信息
• 进行数学运算，包括失量及矩阵操作
• 用if-then-else分支、do循环及用户指令生成执行一系列任务的宏

结构离散化和编号与matlab的一致

APDL常用语法

• !：注释行
• help xx命令：命令行查询xx命令的使用

!%%%%%%%%%%[ANSYS算例]%%%%%begin%%%%
/PREP7 !进入前处理
!=====设置单元和材料
ET,1,SOLID185  !定义1号单元类型为(SOLID185）
MP,EX,1,1e10  !定义1号单元材料弹性模量EX
MP,PRXY,1,0.25  !定义1号单元材料泊松比PRXY

!-----定义8个节点
N,1,0,0,0  !节点1,坐标(0,0,0)，以下类似
N,2,0.2,0,0 
N,3,0,0.8,0 
N,4,0.2,0.8,0
NGEN,2,4,1,4,1,0,0,0.6  !复制节点1~4，新复制的节点编号增量
						!为4，坐标偏移量为（0,0,0.6）,面平移

!------基于节点生成单元
EN,1,1,4,2,6  !由节点1,4,2,6生成单元1，以下类似，共5个单元
EN,2,1,4,3,7 
EN,3,6,7,5,1 
EN,4,6,7,8,4 
EN,5,1,4,6,7 
FINISH


!------进入求解模块
/SOLU
F,7,FZ,-100000,,,  !在节点7处施加FZ，-100000
F,8,FZ,-100000,,,  !在节点8处施加FZ，-100000
D,1,,0,,2,,ALL,,,,,  !对节点1和2，施加固定约束
D,5,,0,,6,,ALL,,,,,  !对节点5和6，施加固定约束
SOLVE !求解
FINISH  !退出该模块


!=====进入一般的后处理模块
/POST1  !进入后处理
PLDISP,1 !计算的变形位移显示(变形前与后的对照)
PRNSOL,U,COMP  !列表给出各个节点位移值
!%%%%%%%%%%[ANSYS算例%%%%% end %%%%

启动软件Mechanical APDL，将上述APDL代码另存为.txt文件后复制，空白处界面右键——粘贴直接导入

所有操作可以在日志文件中查看：list——files——log file

第10讲 有限元方法中的基本性质

节点编号与储存带宽

有限元方法的基本步骤及处理流程：核心是刚度方程

在进行有限元分析时，需要按单元和节点编号所对应的位置，将所形成的单元刚度矩阵装配在整体刚度矩阵中

图示2D连续体第i个单元的节点位移列阵为：
\(\mathbf q^{(i)}=\begin{bmatrix} u_3 & v_3 & u_5 & v_5 & u_6 & v_6\end{bmatrix}^T\)

总共有8个2D节点，即整体刚度矩阵为16*\16

将上述第i个单元的单元刚度矩阵元素填充到整体刚度矩阵中相应位置去，定义第i个单元的半带宽\(d_i\)

• \(d_i\)=(第i个单元中节点编号的最大差值+1)×2 = (6-3+1)×2=8
• 若是3D问题，则\(d_i\)=(第i个单元中节点编号的最大差值+1)×3

整体刚度矩阵的最大半带宽为：\(d_i=max_i\{d_i\};(i=1,2,3,…,n)\)

整体刚度矩阵的存储方式

• 等带宽存储
• 一维变带宽存储

形状函数矩阵与刚度矩阵的性质

以杆单元为例

1. 形状函数的性质

杆单元位移场：
\(u(x)=N_1(x)u_1+N_2(x)u_2=\mathbf N(x) \mathbf q ^e\)

对\(u(x)\)进行考察，特殊情况1

• 左端单位位移，右端固定时
  • 令\(u_1=1,u_2=0\)，有\(u(x)=N_1(x)\)
• 右端单位位移，左端固定时
  • 令\(u_1=0,u_2=1\)，有\(u(x)=N_2(x)\)
• 性质1：\(N_i\)表示在i点的节点位移为1，其它节点位移为0时的单元位移场函数

对\(u(x)\)进行考察，特殊情况2

• 杆单元发生刚体位移\(\bar u_0\)
  • 令\(u(x)=\bar u_0,u_1=u_2=\bar u_0\)。那么\(\bar u_0=N_1(x)\bar u_0+N_2(x)\bar u_0\)
  • 所以，\(1=N_1(x)+N_2(x)\)
• 性质2：单元的形状函数满足\(\sum _{i=1}^n N_i(x) = 1\)，其中n为单元的节点数，它表明形状函数能够描述单元的体位移

2. 单元刚度矩阵系数的性质

单元刚度方程

考虑特殊情况1：主对角线元素

• 左端单位位移，右端固定时
  • 令\(u_1=1,u_2=0\)，有\(k_{11}=p_1\)
• 性质1：单元刚度矩阵的对角线元素\(k_ii\)表示要使单元的第i个节点产生单位位移\(u_i=1\)，而其它节点位移为0时，需在节点i所施加的力

考虑特殊情况2：非主对角线元素

• 右端单位位移，左端固定时

  • 令\(u_1=0,u_2=1\)，有\(k_{12}=p_1\)

• 性质2：单元刚度矩阵的非对角线元素\(k_{ij}(i\ne i)\)表示要使单元的第j个节点产生单位位移\(u_j=1\)，而其它节点位移为0时，需在节点所施加的力

• 性质3：单元刚度矩阵是对称的

  • \(\mathbf{K}^{eT} =[\int_{\Omega^e} \mathbf{B}^T\cdot \mathbf{D} \cdot \mathbf{B}d\Omega]^T=\int_{\Omega^e} \mathbf{B}^T \mathbf{D} \mathbf{B}d\Omega=\mathbf{K}^{e}\)

• 性质4：单元刚度矩阵是半正定的

  • 应变能
  • \(U=\frac{1}{2}\mathbf{q}^{eT}\mathbf{K}^{e}\mathbf{q}^{e}=\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n k_{ij}u_i u_j\)
  • \(\mathbf{q}^{e}\)为非刚体位移时，除非\(\mathbf{q}^{e}=0\)，否则应变能U总为正值
  • \(\mathbf{q}^{e}\)为刚体位移时，有\(\mathbf{q}^{e}\ne 0\)，而此时应变能\(U=0\)，则定有\(\mathbf{K}^{e}=0\)

• 性质5：单元刚度矩阵是奇异的，即\(|\mathbf K^e|=0\)

• 性质6：刚度矩阵的任一行（或列）代表一个平衡力系；当节点位移全部为线位移时（即为\(C_0\)型问题），任一行（或列）的代数和应为零

3. 总体刚度矩阵的性质

• 总体刚度矩阵由单元刚度矩阵组装
  • 对称性
  • 奇异性
  • 半正定性
  • 稀疏矩阵
  • 非零元素显现带状线

边界条件的处理与支反力的计算

整体刚度方程为：\(\mathbf K \mathbf q=\mathbf P\)

两种类型的边界条件

• 零位移：\(\bar{\mathbf q}_a=\mathbf 0\)
• 给定位移：\(\bar{\mathbf q}_a=\bar {\mathbf u}\)

将位移矩阵分为已知的（绿色）和未知的（红色），对应的载荷也是已知的和未知的两部分（已知的位移与未知的载荷是互补关系）

\(\begin{bmatrix} \mathbf K _{aa} & \mathbf K _{ab} \\\mathbf K _{ba} & \mathbf K _{bb}\end{bmatrix}\begin{bmatrix} {\color{Green} \bar {\mathbf q} _{a} } \\ {\color{Red} \bar {\mathbf q} _{b}} \end{bmatrix}=\begin{bmatrix} {\color{Red} \bar {\mathbf P} _{a} } \\ {\color{Green} \bar {\mathbf P} _{b}} \end{bmatrix}\)

1. 直接法

针对零位移的边界条件：\(\bar{\mathbf q}_a=\mathbf 0\)
得到：\(\mathbf K_{bb} \mathbf q_b=\bar{\mathbf P}_b\)
可得：\(\mathbf q_b=\mathbf K_{bb} ^{-1} \bar{\mathbf P}_b\)

针对给定位移的边界条件：\(\bar{\mathbf q}_a=\bar {\mathbf u}\)
得到：\(\mathbf K_{ba} \mathbf q_a+\mathbf K_{bb} \mathbf q_b=\bar{\mathbf P}_b\)
可得：\(\mathbf q_b=\mathbf K_{bb} ^{-1} (\bar{\mathbf P}_b-\mathbf K_{ba}\bar {\mathbf u})\)，继而求出支反力\(\bar{\mathbf P}_a\)

直接法特点：

• a)处理过程直观
• b)待解矩阵规模变小，适合手工处理
• c)矩阵节点编号及排序改变，不利于计算机规范化处理

2. 置一法

主要针对零位移边界条件：\(\bar{ q}_r= 0\)

相应地，对于整体刚度矩阵\(\mathbf K\)的第\(r\)行和第\(r\)列，以及载荷的第\(r\)个元素，使得

• \(k_{rr}=1,k_{rs}=k_{sr}=0(r \ne s)\)
• \(p_r=0\)
  

置一法的特点

• 置1法不应改变原界条件的真实性
  • 由\(k_{rr}\cdot \bar{q}_r=p_r\)可知，\(\bar{ q}_r= 0\)
• 只能处理\(\bar{ q}_r= 0\)的情形
• 保持原矩阵的规模，不需重新排序
• 保持总刚度阵的对称性，利于计算机的规范化处理

3. 乘大数法

针对边界条件：\(\bar{ q}_r= \bar u\)

相应地，对于整体刚度矩阵\(\mathbf K\)的第\(r\)行第\(r\)列的元素\(k_{rr}\)，以及载荷的第\(r\)个元素，使得

• \(k_{rr}\)进行倍乘大数字\(\alpha\)
• \(p_r=\alpha k_{rr} \bar{u}\)

乘大数法的特点

• 置1法不应改变原界条件的真实性
  • 由于\(\alpha\)值很大（如取值\(max|k_{ij}|\times 10^4\)），可得\(\alpha k_{rr}\cdot \bar{q}_r\approx \alpha k_{rr} \bar u\)可知，\(\bar { q}_r \approx \bar u\)
• 既能处理\(\bar{ q}_r= 0\)的情形，又能处理\(\bar{ q}_r= \bar u\)的情形
• 保持原矩阵的规模，不需重新排序
• 保持总刚度阵的对称性，利于计算机的规范化处理

4. 罚函数法

针对边界条件：\({ u}_1= \bar u _1\)
沿着\(u_1\)方向用大刚度\(C\)的弹簧模拟支撑

其应变能为：\(U_s = \frac{1}{2} C(u_1 -\bar u_1)^2\)
系统总势能：\(\Pi =\frac{1}{2}\mathbf{q}^{T}\mathbf{K}\mathbf{q} +\frac{1}{2} C(u_1 -\bar u_1)^2-\mathbf{P}^T\mathbf{q}\)

由\(\partial \Pi /\partial u_i=0(i=1,2,…,n)\)可得刚度方程并展开第一行，有

\((k_{11}+C)u_1+k_{12}u_2+…+k_{1n}u_n=R_1+C\bar{u}_1\)

如果\(C\)值足够大（一般\(C=max|k_{ij}|\times 10^4\)），则有\({ u}_1 \approx \bar u _1\)

继而可求得支反力\(R_1=C(u-\bar u_1)\)

罚函数法的最大好处就是可以直接求出位移边界上的支反力

位移函数构造与收敛性要求

收敛性

• 当节点数或单元插值位移的项数趋于无穷大时，即当单元尺寸趋近于零时，最后的解答如果能够无限地逼近准确解

• 曲线1：收敛于真实解
• 曲线2：收敛于真实解，收敛速度比曲线1慢
• 曲线3：不收敛于真实解
• 曲线4：非单调收敛，不构成真实解的上界或下界
• 曲线5：发散

收敛性的保证

• 单元的内部的收敛：当单元尺寸趋近于零时，使得系统单元的势能至少要存在

  • 位移函数本身应在单元上连续，还要包括至少使得位移函数及对应于应变都为常数的项，即常位移项和常应变项
    

• 单元之间的连接的收敛：位移函数在单元之间应该保证足够的连续，使得“应变”的计算能够存在。

收敛性准则

1. 针对单元内部：常位移项和常应变项

   • 准则1完备性要求：如果在势能泛函中所出现位移函数的最高阶导数是\(m\)阶，则有限元解答收敛的条件之一是单元内的位移场函数的试函数至少是\(m\)阶完全多项式。

2. 针对单元之间：保证足够的连续，使得“应变”的计算能够存在

   • 准则2协调性要求：如果在势能泛函中位移函数出现的最高阶导数是\(m\)阶，则位移试函数在单元交界面上必须具有直至\(m-1\)阶的连续导数，即\(C_{m-1}\)连续性

讨论1：平面问题中，势能函数为

\(\prod =U-W=\frac{1}{2} \int _{\Omega} \sigma _{ij} \epsilon _{ij} d \Omega-[\int _{\Omega} \bar b_i u _i d \Omega+\int _{S_p} \bar p_i u _i dA]\)

\(\epsilon _{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})\)，关于位移的最高阶导数为1，故\(m=1\)

由准则1，形状函数至少应包含完整的一次多项式，即偏导后可以得到常位移项和常应变项

\(\left\{\begin{matrix} u(x,y)=a_0+a_1x+a_2y\\v(x,y)=b_0+b_1x+b_2y\end{matrix}\right.\)

由准则2，位移函数要求\(C_0\)连续，即要求函数本身连续，即仅采用节点位移\((u_i,v_i)\)进行插值后，就可以使得“应变”的计算能够存在。

讨论2：平面梁的弯曲问题中，势能函数为

\(\prod =U-W=\frac{1}{2} \int _{l}EI_z (\frac{d^2v}{dx^2})^2-\int _{l} \bar p(x)\cdot v(x) dx\)

\((\frac{d^2v}{dx^2})^2\)，关于位移的最高阶导数\(m=2\)

由准则1，形状函数至少应包含完整的一次多项式，即偏导后可以得到常位移项和常应变项

\(v(x)=a_0+a_1x+a_2x^2\)

由准则2，位移函数为\(C_1\)连续，即在单元之间的位移函数至少要求一阶导数连续，即必须采用节点位移及节点转角\((v_i,\theta _i)\)进行插值后，才可以使得“应变”的计算能够存在。

因此，基于两个节点的\((v_i,\theta _i)\)，梁单元的插值函数为

\(v(x)=a_0+a_1x+a_2x^2+a_3x^3\)

可以使得\(\epsilon=-y\frac{d^2v}{dx^2}\)的计算存在

常用单元中，能够保证常位移项和常应变项的多项式

• 常轴力杆单元：\(1, x\)
• 平面单元：\(1,x,y\)
• 空间单元：\(1,x,y,z\)
• 平面梁单元：\(1,x, x^2\)
• 平板弯曲单元：\(1, x, y,x^2,xy, y^2\)

单元的位移模式——多项式的选取应考虑：

• （1）唯一确定性：待定系数个数应与节点位移的DOF数相等，选择多项式应从低阶到高阶，尽量选取完全多项式以提高单元的精度
• （2）收敛性准则1(单元内部的完备性）：多项式必须要选择常数项和完备的一次项
• （3）收敛性准则2(单元之间的协调性）：对于\(C_0\)型单元都可以保证，对于\(C_1\)型单元较难保证

可参考由多项式函数构成的 Pasca三角形和上述原则进行函数项次的选取和构造

\(C_0\)单元与\(C_1\)单元

\(C_0\)型单元的续性
\(C_0\)型单元是指势能函中位移函数出现的最高阶导数是1阶，在单元交界面上具有0阶的连续实导数即节点上只要求位移连续。一般的杆单元、平面问题单元、空间问题单元都是\(C_0\)型单元

\(C_1\)型单元的连续性
\(C_1\)型单元是指在势能泛函中位移函数出现的最高阶导数是2阶，在单元交界面上具有1阶的连续导数，即节点上除要求位移连续外，还要求1阶导数连续。梁单元、板单元、壳单元都是\(C_1\)型单元

单元的拼片试验

1. 单元的协调性

• 单元位移函数满足完备性要求→称单元是完备的
• 单元位移函数满足协调性要求→称单元是协调的
• 单元既完备又协调→称单元为协调单元

协调单元的尺寸趋于零时，有限元分析的解答趋于真实解

非协调单元

某些情况下，可以放松对协调性的要求，只要单元能通过拼片试验，有限元分析的解答仍可以收敛于正确解。这种单元称为非协调元

对非协调元，考察其收敛性就是考察其具有常应变的能力。因此，可设计一种数值试验，当对单元片中各节点赋予对应于常应变状态的载荷和位移

节点\(i\)被单元完全包围，节点\(i\)的平衡方程为

\(\sum _{e=1}^m (K^e_{ij}u_j-P^e_i)=0\)

如果能满足，则说明单元能满足常应变要求，当单元尺寸不断减小时，有限元解趋近于真实解

2. 拼片试验

以平面问题中的单元为例，由于对应于常应变的位移是线性位移，即令节点位移为

\(\left\{\begin{matrix} u_j(x,y)=a_0+a_1x_j+a_2y_j\\v_j(x,y)=b_0+b_1x_j+b_2y_j\end{matrix}\right. \qquad (2)\)

若\(i\)点无体积力，无外载荷，则拼片试验的方程为

\(\sum _{e=1}^m K^e_{ij}u_j=0\qquad (3)\)

通过拼片试验的前提是：当赋予各节点以(2)式的位移时，（3)式应该成立

若(3)式不能成立，说明当单元片中各节点具有对应于常应变状态的位移时，节点\(i\)不能保持平衡。则在单元交界面上位移不协调会引起附加应变能

有限元分析数值解的精度与性质

1、求解精度的估计

平面问题为例，单元位移场可展开为

\(\mathbf u=\mathbf u_i +(\frac{\partial \mathbf u}{\partial x})_i\Delta x+(\frac{\partial \mathbf u}{\partial y})_i\Delta y+…\)

• 如果单元尺寸为\(h\)，则上式中的\(\Delta x,\Delta y\)是\(h\)量级
• 如果单元位移函数采用\(p\)阶完全多项式，那么位移解的误差量级\(O(h^{p+1})\)
• 如果应变是位移的\(m\)阶导数给出的，则应变误差量级\(O(h^{(p+1-m)})\)
• 而应变能是应变的平方项表示的，其误差量级为\(O(h^{2(p+1-m)})\)

2、有限元分析结果的下限性质

对象的总势能：\(\Pi =U-W=\frac{1}{2} \mathbf q^T \mathbf K \mathbf q-\mathbf P ^T\mathbf q\)

基于最小势能原理，\(\delta \Pi =0\)，可知刚度方程为\(\mathbf K \mathbf q=\mathbf P\)

平衡时系统的弹性势能，代入刚度方程可知

\(\Pi =\frac{1}{2} \mathbf q^T \mathbf K \mathbf q-\mathbf P ^T\mathbf q=-\frac{1}{2} \mathbf q^T \mathbf K \mathbf q=-U=-\frac{W}{2}\)

最小势能原理的性质：\(\Pi _{approx} \ge \Pi _{exact}\)，即\(U _{approx} \le U_{exact}\)

对应于精确解的刚度矩阵及节点位移

\(\mathbf K _{approx}\mathbf q_{approx}=\mathbf P\)

对应于近似解的刚度矩阵及节点位移

\(\mathbf K _{exact}\mathbf q_{exact}=\mathbf P\)

同一个载荷下，\(-\frac{1}{2} \mathbf q_{approx}^T \mathbf K_{approx} \mathbf q_{approx} \le -\frac{1}{2} \mathbf q_{exact}^T \mathbf K_{exact} \mathbf q_{exact}\)

于是，\(\mathbf q_{approx}^T \mathbf P \le \mathbf q_{exact}^T \mathbf P\)

故，\(|\mathbf q_{approx}| \le |\mathbf q_{exact}|\)

可见，近似解的位移\(\mathbf q_{approx}\)总体上比精确的位移\(\mathbf q_{exact}\)要小，也就是说近似解具有下限性质

有限元模型的性化
有限元方法用有限个自由度近似描述原具有无究穷多个自由度的系统，刚度会增加，系统变得更刚硬，即刚度矩阵总体数值变大。

单元应力计算结果的误差与平均处理

1、应力结果的误差性质

近似解可表达为，精确解+误差，即

\(\left.\begin{matrix} \widehat{\mathbf{u} } =\mathbf{u}_{true}+\delta \mathbf{u}\\\widehat{\mathbf{\epsilon } } =\mathbf{\epsilon }_{true}+\delta \mathbf{\epsilon } \\\widehat{\mathbf{\sigma } } =\mathbf{\sigma}_{true}+\delta \mathbf{\sigma }\end{matrix}\right\}\)

定义应变误差泛函，数学含义：应变差值在加权矩阵作用下的最小二乘

\(\prod (\widehat{\mathbf{\epsilon } },\mathbf{\epsilon }_{true})=\frac{1}{2} \int _{\Omega} (\widehat{\mathbf{\epsilon } }-\mathbf{\epsilon }_{true})^T \cdot \mathbf D \cdot (\widehat{\mathbf{\epsilon } }-\mathbf{\epsilon }_{true}) d \Omega\)
\(=\sum _{e=1}^n \int _{{\Omega}^e} \frac{1}{2} (\widehat{\mathbf{\epsilon } }-\mathbf{\epsilon }_{true})^T \cdot \mathbf D \cdot (\widehat{\mathbf{\epsilon } }-\mathbf{\epsilon }_{true}) d \Omega\)

同理，应变误差泛函

\(\prod (\widehat{\mathbf{\sigma } },\mathbf{\sigma }_{true})=\sum _{e=1}^n \int _{{\Omega}^e} \frac{1}{2} (\widehat{\mathbf{\sigma } }-\mathbf{\sigma }_{true})^T \cdot \mathbf D ^{-1}\cdot (\widehat{\mathbf{\sigma } }-\mathbf{\sigma }_{true}) d \Omega\)

其中，\(\mathbf D\)为加权矩阵

2、Gauss积分点上的应力精度

在GauSs积分点上，应力（应变）的近似解将具有比其他位置高得多的精度。具体对于一个等参元，若采用\(r+1\)个高斯积分点：则高斯积分点上的应变或应力近以解比其它位置的结果具有较高的精度，称该\(r+1\)个高斯积分点是等参元中的最佳应力点。

3、共用节点上应力的平均处理

1.共用节点上应力的直接平均

\(\bar {\sigma}_{kl}(i)=\frac{1}{r}\sum _{e=1} ^r \sigma ^e _{kl} (i)\)

式中\(\sigma ^e _{kl} (i)\)为共用节点让上的平均应力，\(1-r\)为围绕该共用节点周围的全部单元。

2.共用节点上应力的加权平均

由于围绕共用节点周围的客个单元的形状和大小都不一定相同，一种更合理的处理方法是进行加权平均，如果按单元的面积或体积进行加权，则有以下计算公式

\(\bar {\sigma}_{kl}(i)=\sum _{e=1} ^r \eta ^e \sigma ^e _{kl} (i)\)

• 对于2D情形，\(\eta ^e=\frac{A ^e}{\sum _{e=1}^r A ^e}\)
• 对于3D情形，\(\eta ^e=\frac{\Omega ^e}{\sum _{e=1}^r \Omega ^e}\)

以上的处理只是计算结果后处理的一种局部改善，并不能从根本上解决节点应力（应变）精度差的问题

控制误差和提高精度的h方法和p方法

1、h方法：High density

• 不改变各单元上基底函数的配置情况，只通过逐步加密有限元网格来使结果向正确解逼近，此即h方法，也称h-version

该方法可达到一般的工程精度，但由于不采用高阶多项式作基底函数，因而数值稳定性和可靠性都较好，但收敛速度较慢，效率较低

2、p方法：polynomial

• 不改变网格划分情况，只通过增加单元基底函数阶次来使结果向正确解逼近，此即p方法，也称p-version

大量实践表明，p方法的收敛性大大优于h方法。由于p方法使用高阶多项式作基底函数会出现数值稳定性问题

• 多项式的最高阶次：p < 9

3、自适应方法

• 自适应方法运用反馈原理，利用上一步的计算结果来修改有限元模型，其计算量较小，计算精度却得到显著提高

自适应方法选取

• 基于h方法
• 基于p方法
• 基于h-p方法

确定最优网络：对于给定的自由度总数\(N=N_0\)

\(min\{ L(h,p,\lambda)=\Pi _{error}(h,p) -\lambda (\int _{\Omega} n(h,p)d\Omega-N_0)\}\)

• \(h\)：单元特质尺寸
• \(p\)：多项式阶次
• \(\lambda (\int _{\Omega} n(h,p)d\Omega-N_0)\)：采用Lagrange乘子法

自适应方法流程：

1. 事前误差估计
2. 初始数值模拟
3. 误差分析及估计
4. 选择 h,p
5. 改进方案：采用高精度的数值分析方法
6. 方程的求解
7. 事后误差估计，是否满足精度要求
   1. No，回到步骤3
   2. Yes，继续执行8
8. 获得满意的结果

第11讲 高阶及复杂单元

1D 高阶单元

1、自然坐标（1D）

定义，\(L_1=l_1/l,L_2=l_2/l,L_1+L_2=1\)，P点坐标为\((L_1,L_2)\)

2、一次杆单元（基于自然坐标）

节点位移列阵：\(\mathbf q^e=[u_1,u_2]^T\)
单元位移模式：
\(u(x)=a_1+a_2x=a_1'L_1+a_2'L_2=N_1u_1+N_2u_2=\mathbf N(L_1,L_2)\mathbf q^e\)

形状函数表达(x坐标)：\(\mathbf N=[N_1,N_2],N_1=1-x/l,N_2=x/l\)
形状函数表达(自然坐标)：\(\mathbf N=[N_1,N_2],N_1=L_1,N_2=L_2\)

3、二次杆单元

3个节点

节点位移列阵：\(\mathbf q^e=[u_1,u_2,u_3]^T\)
单元位移模式：
\(u(x)=a_1+a_2x+a_3x^2=a_1'L_1+a_2'L_2+a_3'L_1L_2\)
\(=N_1u_1+N_2u_2+N_3u_3=\mathbf N(L_1,L_2)\mathbf q^e\)

形状函数表达(x坐标)：\(\mathbf N=[N_1,N_2,N_3],\)
\(N_1=(1-2\frac{x}{l})(1-\frac{x}{l}),N_2=4\frac{x}{l}(1-2\frac{x}{l}),N_3=-\frac{x}{l}(1-2\frac{x}{l})\)
形状函数表达(自然坐标)：\(\mathbf N=[N_1,N_2,N_3],\)
\(N_1=L_1(2L_1-1),N_2=4L_1L_2,N_3=L_2(2L_2-1)\)

4、三次杆单元

4个节点

节点位移列阵：\(\mathbf q^e=[u_1,u_2,u_3,u_4]^T\)
单元位移模式：
\(u(x)=a_1+a_2x+a_3x^2+a_4x^4=a_1'L_1+a_2'L_2+a_3'L_1L_2+a_4'L_1^2L_2\)
\(=N_1u_1+N_2u_2+N_3u_3+N_4u_4=\mathbf N(L_1,L_2)\mathbf q^e\)

形状函数表达(x坐标)：\(\mathbf N=[N_1,N_2,N_3,N_4],\)
\(N_1=(1-\frac{3x}{l})(1-\frac{3x}{2l})(1-\frac{x}{l})\)
\(N_2=9\frac{x}{l}(1-\frac{3x}{2l})(1-\frac{x}{l})\)
\(N_3=-\frac{9x}{2l}(1-\frac{3x}{l})(1-\frac{x}{l})\)
\(N_4=\frac{x}{l}(1-\frac{3x}{l})(1-\frac{3x}{2l})\)
形状函数表达(自然坐标)：\(\mathbf N=[N_1,N_2,N_3,N_4],\)
\(N_1=L_1(1-\frac{9}{2}L_1L_2)\)
\(N_2=-\frac{9}{2}L_1L_2(1-3L_1)\)
\(N_3=9L_1L_2(1-\frac{3}{2}L_1)\)
\(N_4=L_2(1-\frac{9}{2}L_1L_2)\)

5、一般高阶杆单元

具有n个节点的1D杆单元

\(u(x)=\sum _{i=1}^n N_iu_i\)

根据形状函数的性质：

• \(N_i(x_j)=\delta _{ij}\)
• \(\sum _{i=1}^n N_i=1\)
• 使用Lagrange插值公式，\(N_i=l_i^{(n-1)}(\xi)=\Pi _{j=1,j\ne i}^n \frac{\xi -\xi _j}{\xi _i- \xi _j},\xi=\frac{x-x_i}{l}\)

6、Hermite单元(一般高阶\(C_1\)型单元)

Hermite多项式进行函数插值，2节点梁单元（要求\(C_1\)连续）

2D 高阶单元

1、自然坐标（2D）

定义面积坐标，\(L_i=A_i/A,L_j=A_j/A,L_m=A_m/A,L_i+L_j+L_m=1\)，P点坐标为\((L_i,L_j,L_m)\)

2、平面3节点三角形单元（基于自然坐标）

节点位移列阵：\(\mathbf q^e=[u_1,v_1,u_2,v_2,u_3,v_3]^T\)
单元位移模式：
\(u(x,y)=\bar a_0+\bar a_1x+\bar a_2y\)
\(=\bar a_0'L_1+\bar a_1'L_2+\bar a_2'L_3\)
\(=N_1u_1+N_2u_2+N_3u_3\)
\(=\mathbf N(L_1,L_2,L_3)\mathbf q^e\)

形状函数表达(自然坐标)：\(\mathbf N=[N_1,N_2,N_3],N_1=L_1,N_2=L_2,N_3=L_3\)

\(v(x,y)\)位移模式同理

3、6节点三角形二次单元

节点位移列阵：\(\mathbf q^e=[u_1,v_1,…,u_6,v_6]^T\)
单元位移模式：
\(u(x,y)=a_1+a_2x+a_3y+a_4x^2+a_5xy+a_6y^2\)
\(= a_1'L_1+a_2'L_2+a_3'L_3+a_4'L_1L_2+a_5'L_2L_3+a_6'L_3L_1\)
\(=N_1u_1+N_2u_2+N_3u_3+…+N_6u_6\)
\(=\mathbf N(L_1,L_2,…，N_6)\mathbf q^e\)

形状函数表达（自然坐标）
\(N_1=L_1(2L_1-1)\)
\(N_2=L_2(2L_2-1)\)
\(N_3=L_3(2L_3-1)\)
\(N_4=4L_1L_2\)
\(N_5=4L_2L_3\)
\(N_6=4L_3L_1\)

\(v(x,y)\)位移模式同理

4、10节点三角形三次单元

节点位移列阵：\(\mathbf q^e=[u_1,v_1,…,u_{10},v_{10}]^T\)
单元位移模式：
\(u(x,y)=a_1+a_2x+a_3y+a_4x^2+a_5xy+a_6y^2+a_7x^3+a_8x^2y+a_9xy^2+a_{10}y^3\)
\(= a_1'L_1+a_2'L_2+a_3'L_3+a_4'L_1L_2+a_5'L_2L_3+a_6'L_3L_1+a_7L_1^2L_2+a_8'L_2^2L_3+a_9'L_3^2L_1+a_{10}'L_1L_2L_3\)
\(=N_1u_1+N_2u_2+N_3u_3+…+N_{10}u_{10}\)
\(=\mathbf N(L_1,L_2,…，N_{10})\mathbf q^e\)

形状函数表达（自然坐标）
\(N_i=\frac{L_i-\frac{2}{3}}{\frac{1}{3}}\cdot \frac{L_i-\frac{1}{3}}{\frac{2}{3}} \cdot \frac{L_i}{1},(i=1,2,3)\)
\(N_{10}=3L_1\cdot 3L_2\cdot 3L_3=27L_1L_2L_3\)

\(N_4=\frac{9}{2}L_1L_2(3L_1-1)\)
\(N_5=\frac{9}{2}L_1L_2(3L_2-1)\)
\(N_6=\frac{9}{2}L_2L_3(3L_2-1)\)
\(N_7=\frac{9}{2}L_2L_3(3L_3-1)\)
\(N_8=\frac{9}{2}L_1L_3(3L_3-1)\)
\(N_9=\frac{9}{2}L_1L_3(3L_1-1)\)

\(v(x,y)\)位移模式同理

5、基于Lagrange插值的矩形单元

具有(r+1)列和(p+1)行节点的矩形单元

单元位移模式

\(u(\xi,\eta)=a_1+a_2\xi+a_3\eta+…+a_i\xi ^r\eta ^p\)

在\(\xi\)方向的\((r+1)\)列节点中，构造一个插值函数在第\(I\)个节点上等于1，而在其他节点上等于0，则该函数为

\(l^{(r)}_I(\xi)=\Pi _{j=0}^r \frac{\xi -\xi {j}}{\xi _I -\xi _j}\)

同理在\(\eta\)方向也可以得到插值函数

\(l^{(p)}_J(\eta)=\Pi _{j=0}^p \frac{\eta -\eta {j}}{\eta _J -\eta _j}\)

对以上两个方向的Lagrange多项式进行乘积运算可得到节点\(i(I,J)\)的插值函数\(N_i\)为

\(N_i=l^{(r)}_I(\xi)l^{(p)}_J(\eta)\)，满足\(N_i(x_j)=\delta _{ij}\)

该单元在每一边界上的节点数和插值函数在边界上的变化是协调的这也保证了单元之间函数的协调性

对于线性、二次和三次函数变化的Lagrange单元

• 随着函数插值阶次的增高，必然要增加内部节点，但这些节点自由度的增加一般不能显著提高单元的精度

Lagrange矩形单元中的多项式

• 完全多项式：利于提高单元的精度
• 非完全多项式：对于提高单元的精度作用不大

主要取边节点的Serendipit单元在实际应用中得到比Lagrange单元更广泛

6、Serendipity矩形单元：尽可能在边界上取节点的高阶单元

第1步 假定只有4个角节点此时的形状函数为

\(N_i'=\frac{1}{4}(1+\xi _i\xi)(1+\eta _i \eta),(i=1,2,3,4)\)

第2步 在1-2边中点增加节点5构造节点5的形状函数为

\(N_5=\frac{1}{2}(1-\xi ^2)(1- \eta)\)

此时的位移场函数为

\(u(x,y)=N_1'u_1+N_2'u_2+N_3'u_3+N_4'u_4+N_{5}u_{5}\)

第3步增加其他边的内节点6，7，8，进行类似补偿计算，可得

\(N_1=N_1'-\frac{1}{2}N_5-\frac{1}{2}N_8\)
\(N_2=N_2'-\frac{1}{2}N_5-\frac{1}{2}N_6\)
\(N_3=N_3'-\frac{1}{2}N_6-\frac{1}{2}N_7\)
\(N_4=N_4'-\frac{1}{2}N_7-\frac{1}{2}N_8\)

\(N_5=\frac{1}{2}(1-\xi ^2)(1- \eta)\)
\(N_6=\frac{1}{2}(1+\xi)(1- \eta ^2)\)
\(N_7=\frac{1}{2}(1-\xi ^2)(1+ \eta)\)
\(N_5=\frac{1}{2}(1-\xi )(1- \eta ^2)\)

这些函数满足， \(N_i(P_j)=\delta _{ij}\)， \(\sum _{i=1}^n N_i=1\)

6、Serendipity矩形单元

三次Serendipity单元与三次Lagrange矩形单元的函数项次比较

Serendipity单元中对于完全多项式以外的高次项使用得较少。这使得当节点数一定时 Serendipity单元的精度较高。

3D 高阶单元

1、10节点四面体二次单元

节点位移列阵：\(\mathbf q^e=[u_1,v_1,w_1,…,u_{10},v_{10},w_{10}]^T\)
单元位移模式：对3D问题，由Pascal三角形可知，一个完备的二次函数有10项。由此，构造具有10个节点的四面体，包含4个角节点、6个分布在6条棱边上的中点

\(u(x,y,z)=a_1+a_2x+a_3y+a_4z+a_5xy\)
\(+a_6yz+a_7xz+a_8x^2+a_9y^2+a_{10}z^2\)

基于自然坐标（体积坐标）的形状函数

\(N_i=(2{L_i}-1)L_i,(i=1,2,3,4)\)
\(N_{5}=4L_1L_2\)
\(N_{6}=4L_2L_3\)
\(N_{7}=4L_1L_3\)
\(N_{8}=4L_1L_4\)
\(N_{9}=4L_2L_4\)
\(N_{10}=4L_3L_4\)

2、20节点四面体三次单元

节点位移列阵：\(\mathbf q^e=[u_1,v_1,w_1,…,u_{20},v_{20},w_{20}]^T\)
单元位移模式：对3D问题，由Pasca三角形可知，一个完备的三次函数有20项。由此，构造具有20个节点的四面体，包含4个角节点、12个分布在6条棱边上的三等分节点、4个面心节点

\(u(x,y,z)=a_1+a_2x+a_3y+a_4z+a_5xy\)
\(+a_6yz+a_7xz+a_8x^2+a_9y^2+a_{10}z^2\)
\(+a_{11}x^3+a_{12}y^3+a_{13}z^3+a_{14}x^2y+a_{15}xy^2+a_{16}y^2z\)
\(+a_{17}yzx^2+a_{18}x^2z+a_{19}xz^2+a_{20}xyz\)

基于自然坐标（体积坐标）的形状函数

角节点
\(N_i=\frac{1}{2}L_i(3{L_i}-1)(3{L_i}-2),(i=1,2,3,4)\)

棱边上的三等分点
\(N_{5}=\frac{9}{2}L_1L_2(3L_1-1)\)
\(N_{6}=\frac{9}{2}L_1L_2(3L_2-1)\)
\(N_{7}=\frac{9}{2}L_2L_3(3L_2-1)\)
\(N_{8}=\frac{9}{2}L_2L_3(3L_3-1)\)
\(N_{9}=\frac{9}{2}L_1L_3(3L_3-1)\)
\(N_{10}=\frac{9}{2}L_1L_3(3L_1-1)\)
……

面心节点
\(N_{17}=27L_1L_2L_4\)
\(N_{18}=27L_2L_3L_4\)
\(N_{19}=27L_1L_3L_4\)
\(N_{20}=27L_1L_2L_3\)

3、Lagrange正六面体高阶单元

与构造2D问题Lagrange矩形单元的插值函数类似，该单元的插值函数直接由三入坐标方向的Lagrange插值多项式的乘积来获得。

单元位移模式

\(u(\xi,\eta,\zeta)=a_1+a_2\xi+a_3\eta+a_4\zeta+…+a_i\xi ^m\eta ^n \zeta ^p\)

单元的形状函数

\(N_i=l^{(m)}_I(\xi)l^{(n)}_J(\eta) l^{(p)}_K(\zeta)\)，满足\(N_i(x_j)=\delta _{ij}\)

其中\(m,n,p\)分别代表每一坐标方向的节点划分数减1，也即为每一坐标方向Lagrange多项式的次数\(I,J,K\)表示节点\(i\)在每一坐标方向的行列式号。

4、Serendipity正六面体单元

与构造Serendipity矩形单元的形状函数类似，同样可以构造出各种节点的 Serendipity正六面体单元

基于薄板理论的弯曲板单元

1、基本变量与方程

• 板壳结构中的几特点是其厚度远小于其他两方向的尺寸
• 引入一定假设对厚度方向的受力特点进行简化，这些简化假设叫做Kirchhoff假定
• 由于薄板中要保持转角的连续，可以承受弯矩，因而薄板问题是\(C_1\)问题

Kirchhoff假定

• 薄板中面法线变形后仍为直法线，具厚度方向正应变很小，有
  • \(\gamma _{zx}=0,\gamma _{zy}=0,\epsilon _{z}=0\)
• 假设薄板中面无横向位移，则会导出
  • \(u(x,y,z=0)=0,v(x,y,z=0)=0\)
  • 进一步得到，\(u(x,y,z)=-z\frac{\partial w}{\partial x},v(x,y,z)=-z\frac{\partial w}{\partial y}\)
• 应力\(\sigma _{zz}\)引起的形变很小，可以忽略

三大类变量

位移：\(w(x,y)=w(x,y,z=0),u(x,y,z)=-z\frac{\partial w}{\partial x},v(x,y,z)=-z\frac{\partial w}{\partial y}\)
应变：\(\epsilon_{x}=-z\frac{\partial ^2 w}{\partial x^2},\epsilon_{y}=-z\frac{\partial ^2 w}{\partial y^2},\gamma_{xy}=-2z\frac{\partial ^2 w}{\partial xy}\)
应力：\(\sigma _x=\frac{12M_y}{h^3}z,\sigma _y=\frac{12M_x}{h^3}z,\tau _{xy}=\frac{12M_{xy}}{h^3}z\)

三大类方程

• 平衡方程

\(D_0(\frac {\partial ^4 w }{\partial x^4}+2\frac {\partial ^4 w }{\partial x^2y^2}+\frac {\partial ^4 w }{\partial y^4})=\bar p(x,y)\)
或者，\(\frac {\partial ^4 M_x }{\partial x^2}+2\frac {\partial ^2 M_{xy} }{\partial xy}+\frac {\partial ^2 M_y }{\partial y^2})+\bar p(x,y)=0\)
其中，\(D_0=\frac{Eh^3}{12(1-\mu ^2)}\)为板弯曲刚度

• 物理方程

\(\left\{\begin{matrix} \sigma _x=\frac{E}{1-\mu ^2}(\epsilon _x+\mu \epsilon _y) \\ \sigma _y=\frac{E}{1-\mu ^2}(\epsilon _y+\mu \epsilon _x)\\ \tau _{xy}=\frac{E}{2(1+\mu )}\gamma _{xy})\end{matrix}\right.\)

或者，\(\mathbf{M}=\widetilde{\mathbf{D}} \mathbf {\kappa }\)
其中，广义力\(\mathbf{M}=\begin{bmatrix} M_x & M_y & M_{xy}\end{bmatrix}^T\)

位移\(\mathbf {\kappa }=\begin{bmatrix} \epsilon_x & \epsilon_y & \gamma_{xy}\end{bmatrix}^T\)
弹性系数矩阵\(\widetilde{\mathbf{D}}=\frac{E}{12(1-\mu^2)}\begin{bmatrix} 1 & \mu & 0\\ \mu & 1 & 0\\ 0 & 0 & \frac{1-\mu}{2} \end{bmatrix} =D_0 \begin{bmatrix} 1 & \mu & 0\\ \mu & 1 & 0\\ 0 & 0 & \frac{1-\mu}{2} \end{bmatrix}\)

• 几何方程

参见应变表达式\(\epsilon_{x}=-z\frac{\partial ^2 w}{\partial x^2},\epsilon_{y}=-z\frac{\partial ^2 w}{\partial y^2},\gamma_{xy}=-2z\frac{\partial ^2 w}{\partial xy}\)

边界条件

• 位移和转角边界条件

\(\left.\begin{matrix}w=\bar w \\\frac{\partial w}{\partial n} =\bar \theta\end{matrix}\right\}\quad on \quad S_1\)

其中，\(\bar w,\bar \theta\)为给定的挠度和转角，\(n\)为边界法线方向

• 位移和力矩边界条件

\(\left.\begin{matrix}w=\bar w \\ M _n =\bar M _n\end{matrix}\right\}\quad on \quad S_2\)

其中，\(\bar w,\bar M _n\)为给定的挠度和力矩，\(n\)为边界法线方向

• 位移和集中剪力边界条件

\(\left.\begin{matrix} M _n =\bar M _n \\ Q_n + \frac {\partial M_{nm}}{\partial m} =\bar Q_n \end{matrix}\right\}\quad on \quad S_3\)

其中，\(\bar M_n,\bar Q _n\)为给定的力矩和横向载荷，\(n\)为边界法线方向，\(m\)为边界切向，\(Q_n\)为边界界面上单位长度的横向剪力

总势能为
\(\Pi =\int (\frac{1}{2} \mathbf{\kappa }^T\mathbf{\widetilde{D}}\mathbf{\kappa }-\bar{p}w)dxdy-\int _{S_3}\bar{Q}_nw ds-\int _{S_1+S_2}\bar{M}_n \frac{\partial w}{\partial n} ds\)

2、4节点矩形非协调薄板单元

节点位移列阵，\(\mathbf {q}^e=\begin{bmatrix} \mathbf {q}_1^T & \mathbf {q}_2^T &\mathbf {q}_3^T &\mathbf {q}_4^T \end{bmatrix}^T\)

每个节点的自由度
\(\mathbf {q}_i= \begin{bmatrix} w_i &\frac{ \partial w _i}{\partial y} & -\frac{ \partial w _i}{\partial x}\end{bmatrix}^T=\begin{bmatrix} w_i & \theta _{xi} & \theta _{yi} \end{bmatrix}^T,(i=1,2,3,4)\)

其中，\(\theta _{xi},\theta _{yi}\)分别为绕\(x\)轴和\(y\)轴的转角

单元位移模式，使用多项式插值（使之满足完备性要求）

\(w(x,y)=a_1 +a_2x+a_3y+a_4x^2 +a_5xy+a_6y^2\)
\(+a_7x^3 + a_8x^2y+ a_9xy^2 + a_{10}y^3 +a_{11}x^3y+a_{12}xy^3\)

单元交界面上\(w\)是连续的，但单元之间的法向导数连续性一般不能满足，所以这种单元是非协调的。但由于这种单元能通过拼片实验，所以当单元网格不断缩小时，计算结果还是可以收敛于精确解的。

3、3节点三角形非协调薄板单元

4、常用弯曲薄板单元一览

4节点矩形薄板单元

3节点三角形薄板单元

子结构与超级单元

结构的重复性

• 几何空间上：采用子结构
• 计算时间上：采用超级单元

1、子结构：系统中具有相同特征和性质的局部结构

子结构划分及应用的方法

• （1）取具有重复性的结构作为子结构（可以有多级子结构）

• （2）对最底层子结构进行分析，形成刚度方程并缩聚

  设第\(k\)级子结构刚度方程为

  

  进行缩聚处理
  

• （3）将子结构进行拼装形成上一级子结构

• （4）对多级子结构全部处理后，得最终的整体刚度方程，求

• （5）将结果回代，再求各级子结构内部的节点位移和其它物理量

2、超级单元

主要用于多次的迭代计算（时间历程）

• 超级单元：一种广义的特定单元，是一个经缩聚内部节点自由度后的子结构，它只有与外部有连接关系的节点位移自由度
• 目的：减小计算量，特别是在需要多次选代的复杂计算中有明显优越性
• 效果：大大减小每次生成刚度矩阵的计算量，同时也减小了计算规模，获得较高的计算效率

使用刚度方程描述

对从节点的节点位移\(\mathbf{q}_s\)进行缩聚

方程（9）就代表该超级单元的单元刚度方程，\(\widetilde{\mathbf K}\)就是该超级单元的刚度矩阵， \(\widetilde{\mathbf Q}\)为超级单元的外载节点列阵

第12讲 有限元分析的应用领域引论(1)

结构振动的有限元分析：基本原理

1、结构振动问题的概述

任何变形体都存在的固有频率和振动模态，当有外界的激振作用时，会产生一系列的响应，除结构的静力分析外，结构的振动分析也是结构评价的另一个重要方面，对结构的工作状态及功能控制具有重要意义。

2、结构振动问题的基本变量

三大变量（均为坐标位置\(\xi (x,y,z)\)和时间\(t\)的函数）

2D问题

• 位移\(u_i(\xi,t) \to\)位移\(u(\xi,t),v(\xi,t)\)
• 应变\(\epsilon_{ij}(\xi,t) \to\)应变\(\epsilon_x(\xi,t),\epsilon_y(\xi,t),\gamma_{xy}(\xi,t)\)
• 应力\(\sigma_{ij}(\xi,t) \to\)应力\(\sigma_x(\xi,t),\sigma_y(\xi,t),\tau_{xy}(\xi,t)\)

3、结构振动问题的基本方程

（1）平衡方程

微小体元dxdydz在动力学状态下的平衡关系，由D' Alembert原理

\(\sigma_{ij,j}(t)-\bar b_{i}(t)-\rho \ddot{u}_i(t) -v\dot{u}_i(t)=0\)

其中\(\rho\)为密度，\(v\)为阻尼系数，\(\bar b_i(t)\)为所作用的体积力，\(\ddot{u}_i(t),\dot{u}_i(t)\)分别为\(u_i(t)\)对时间\(t\)的二次导数和一次导数，即表示\(i\)方向的加速度和速度

（2）几何方程

\(\epsilon _{ij}(t)=\frac{1}{2}(u_{i,j}(t)+u_{j,i}(t))\)

（3）物理方程

\(\sigma _{ij}(t)=D_{ijkl}\epsilon _{kl}(t)\)

（3）边界条件

\(u_i(t)=\bar u_{i}(t) \quad On \quad Su\)，位移边界条件\(BC(u)\)
\(\sigma_{ij}(t)n_j=\bar p_{i}(t) \quad On \quad Sp\)，力边界条件\(BC(p)\)

特别的，初始条件IC(initial condition）
\({u}_i (\xi,t =0) =\bar u_i(\xi)\)
\(\dot u _i(\xi,t =0)=\dot {\bar u}_i(\xi)\)

4、结构振动问题的虚功原理

在引入惯性力和阻尼力的基础上，写出平衡方程及力边界条件的等效积分形式

\(\delta \Pi =\int _\Omega -[\sigma _{ij,j}-\rho \ddot{u}_i -v\dot{u}_i+\bar b_i]\delta u_id\Omega+\int _{S_p}[\sigma _{ij}n_j-\bar p_i]\delta u_idA=0\)

对上述方程右端的第一项进行分部积分(应用Gauss-Green公式），有

这就是动力学问题的虚位移方程（虚功方程）

5、结构振动问题的有限元分析列式

单元的节点位移矩阵：\(\mathbf q_t^e(t)=\begin{bmatrix} u_1(t) & v_1(t) & w_1(t) & … & u_n(t) & v_n(t) & w_n(t)\end{bmatrix}^T\)

单元内的位移插值函数：\(\mathbf u^e(\xi,t)=\mathbf N(\xi)\cdot \mathbf q^e_t(t)\)

应变、应力、速度、加速度由节点位移矩阵表示为

\(\mathbf {\epsilon}^e(\xi,t)=[\partial]\mathbf {u}^e=[\partial]\mathbf N(\xi) \mathbf q^e_t(t)=\mathbf B(\xi)\cdot \mathbf q^e_t(t)\)
\(\mathbf {\sigma}^e(\xi,t)=\mathbf D \cdot \mathbf {\epsilon}^e=\mathbf D \cdot \mathbf B(\xi)\cdot \mathbf q^e_t(t)=\mathbf S(\xi)\cdot \mathbf q^e_t(t)\)
\(\dot {\mathbf u}^e(\xi,t)=\mathbf N(\xi)\cdot \dot{\mathbf q}^e_t(t)\)
\(\ddot {\mathbf u}^e(\xi,t)=\mathbf N(\xi)\cdot \ddot{\mathbf q}^e_t(t)\)

将相关物理量代入前面的虚功方程

节点位移微元\(\delta {\mathbf q}_t^e(t)\)具有任意性，故

\(\mathbf {M}^e \ddot {\mathbf q}_t^e(t)+\mathbf {C}^e \dot {\mathbf q}_t^e(t)+\mathbf {K}^e {\mathbf q}_t^e(t)=\mathbf {P}^e_t(t)\)

其中，
单元的质量矩阵\({\mathbf M}^e=\int _{\Omega _e} \rho {\mathbf N} ^T{\mathbf N}d\Omega\)
单元的阻尼矩阵\({\mathbf C}^e=\int _{\Omega _e} v {\mathbf N} ^T{\mathbf N}d\Omega\)
单元的刚度矩阵\({\mathbf K}^e=\int _{\Omega _e} {\mathbf B} ^T{\mathbf D}{\mathbf B}d\Omega\)
\({\mathbf P}^e_t=\int _{\Omega _e} {\mathbf N} ^T \bar{\mathbf b}d\Omega+\int _{S _p} {\mathbf N} ^T \bar{\mathbf p}dA\)

进行单元的装配，得到结构振动总刚度方程

\(\mathbf {M} \ddot {\mathbf q}_t+\mathbf {C} \dot {\mathbf q}_t+\mathbf {K} {\mathbf q}_t=\mathbf {P}_t\)

6、结构振动问题的讨论

（1）静力学情形(static case）与时间无关，退化为\(\mathbf {K} {\mathbf q}=\mathbf {P}\)

（2）无阻尼情形（undamped case）,\(v=0\)，退化为\(\mathbf {M} \ddot {\mathbf q}_t+\mathbf {K} {\mathbf q}_t=\mathbf {P}_t\)

（3）无阻尼自由振动(free vibration of undamped system），自由振动
\(v=0,\mathbf {P}_t=0\)，退化为\(\mathbf {M} \ddot {\mathbf q}_t+\mathbf {K} {\mathbf q}_t=0\)
解的形式为简谐振动，\(w\)为常数
\({\mathbf q}_t=\widetilde{\mathbf q}_t \cdot e^{iwt}\)
有非零解，特征方程\(|\mathbf{K}-\omega ^2\mathbf{M}|=0\)

• 特征值\(\lambda =\omega ^2\)
• 自然圆频率(\(rad/s\))\(\omega\)
• 自然频率\(f_{rq}=\frac{\omega }{2 \pi}\)
• 特征向量\(\widetilde{\mathbf q}\)对应振动频率\(\omega\)的振型

7、结构振动问题中的质量矩阵

\(\mathbf {M}^e \ddot {\mathbf q}_t^e(t)+\mathbf {C}^e \dot {\mathbf q}_t^e(t)+\mathbf {K}^e {\mathbf q}_t^e(t)=\mathbf {P}^e_t(t)\)

• 一致质量矩阵
  直接由形状函数矩阵推导质量矩阵，“一致”：这里用的形状函数与推导刚度矩阵所用的形状函数是一致的，其中各系数之间存在耦合
• 集中质量矩阵
  将质量集中到节点上，“集中”：系数都集中在矩阵对角线上，各个系数相互独立、解耦，便于求解

杆单元的质量矩阵

梁单元的质量矩阵

质量矩阵类型	梁单元的质量矩阵	平面三节点三角形单元	平面四节点四角形单元

结构振动的有限元分析实例

一个阶梯杆结构的轴向自由振动分析

进行轴向自由振动频率和模态分析

• 弹性模量：\(E\)
• 密度：\(\rho\)
• 横截面积关系：\(A^{(1)}=2A^{(2)}=2A\)

1、单元离散

离散方式：自然离散

单元编号

代入总体的自由振动方程（\(w\)为自然圆频率）

自由振动的特征方程

\(|\mathbf{K}-\omega ^2\mathbf{M}|=0\)

令\(\lambda ^2=\frac{\rho L^2 \omega ^2}{24E}\)，得到

\(18\lambda ^2(1-2\lambda ^2)(\lambda ^2-2)=0\)

振型\(\widehat{q} _1,\widehat{q} _2,\widehat{q} _3\)关于质量矩阵\(M\)正交，即

\(\widehat{\mathbf q} _i^T \mathbf M \widehat{\mathbf q} _j=\left\{\begin{matrix}a_i \quad 当i=j \\0 \quad 当i\ne j\end{matrix}\right.\)

弹塑性问题的有限元分析：基本原理

1、弹塑性问题的概述
弹塑性问题(elastic-plasticity problem)是指变形体的力学行为呈现出超出弹性极限后的塑性行为。

在大多数情况下，人们在进行结构设计时一般都以材料的弹性极限作为依据，还考虑有一定的安全系数，也就是说，处于工作状态中的材料应该处于弹性范围。

2、弹塑性问题的物理方程

研究弹塑性问题(elastic-plastic problem)的关键在于物理方程的处理
材料弹塑性行为实验（单向拉伸或压缩）的特征

1. 初始屈服条件，即屈服准则
2. 材料屈服后，塑性应变增长的规律，即流动法则
3. 新的屈服极限，即强化准则

（1）屈服准则——确定材料产生屈服时的临界应力状态(3D）
大量的实验表明：多数材料的塑性屈服(plastic yielding)与静水压力无关
可以利用等倾八面体

定义等效应力

初始屈服条件

定义屈服函数

（2）塑性流动准则——确定塑性应变分量在塑性变化时的大小、方向

塑性状态的后继行为判断

（3）强化准则——描述屈服面如改变，确定新的屈服面状态

• 等向强化(isotropic hardening)模型
• 随动强化(kinematic hardening)模型
• 混合强化(非等向)(anisotropic hardening)模型

在发生塑性强化的情况下，材料的临界屈服应力将随看塑性应变的积累而发生变化

材料塑性行为的几种典型的模型

• a.双线性随动强化（bilinear kinematic）
• b.多线性随动强化（multilinear kinematic）
• c.双线性等向强化（bilinear isotropic）
• d.多线性等向强化(multilinear isotropic)
• e.非等向强化(Anisotropic)
• f.Drucker-Prager模型

• 基于全量理论的有限元分析列式
  假设：整个加载过程为比例加载(proportionally loading），其结果只与状态有关，与加载过程无关

• 基于增量理论的有限元分析列式
  增量理论考虑真实的加载过程，即变形结果与加载历史有关

弹塑性问题的有限元分析：非线性问题求解

1、非线性方程求解的Newton-Raphson（N-R）选代法

有限元列式——非线性方程组：\(\mathbf K ^{ep}(\mathbf q)\cdot \Delta \mathbf q=\Delta \mathbf P\)

非线性方程组求解

• 直接选代法
• Newton-Raphson（N-R）迭代法
• 改进的N-R迭代法

Newton-Raphson（N-R）迭代法

• 思想：分步逼近计算
• 主要步骤：将总载荷分成一系列载荷段，在每一载荷段内进行非线性方程的迭代
• 特点：在迭代中非线性方程变为线性方程

1.将总载荷分为一系列载荷段

\(\mathbf P=\mathbf P^{(1)}+\mathbf P^{(2)}+\mathbf P^{(3)}+…+\mathbf P^{(n)}\)

2. 多步迭代

在载荷段\((\mathbf P^{(k)},\mathbf P^{(k+1)})\)内保证收敛，再进入下一个载荷段

第\(k\)个载荷步，第\(i\)次迭代

\(\mathbf K ^{ep}_T(\mathbf q_i^{(k)})\cdot \Delta \mathbf q _i ^{(k)}=\Delta \mathbf P_i ^{(k)}\)

其中，
\(\mathbf K ^{ep}_T(\mathbf q_i^{(k)})\)为弹塑性切线的刚度矩阵，对应\(\mathbf K ^{ep}_{\tau}(\mathbf q)\)
\(\Delta \mathbf P_i ^{(k)}=\mathbf P^{(k)}_{i+1}-\mathbf P^{(k)}_i\)

解出\(\Delta \mathbf q _i ^{(k)}\)

载荷段内迭代到收敛

3. 所有载荷段：循环迭代，累加结果

2、修正的Newton-Raphson（N-R）迭代法

• 常规的N-R迭代：需要每次重新形成切线刚度矩阵并求逆，计算量很大
• 改进的N-R迭代：始终采用初始的切线刚度矩阵，保持不变大大减少计算量

第13讲 有限元分析的应用领域引论(2)

传热问题的有限元分析：基本原理

1、传热问题的概述
传热(heat transfer)是日常生活和工程实际中广泛存在的自然现象

进行热题的分析一般包括两部分内容：

• 1、通过传热分析来确定温度场
• 2、在获得温度场的基础上，计算所产生的热应力

2、传热问题的控制方程

• Fourier传热定律
• 能量守恒定律

物体的瞬态温度场\(T(x,y,z,t)\)

\(\frac{\partial }{\partial x} (\kappa _x \frac{\partial T}{\partial x} )+\frac{\partial }{\partial y} (\kappa _y \frac{\partial T}{\partial y} )+\frac{\partial }{\partial z} (\kappa _z \frac{\partial T}{\partial z} )+\rho Q=\rho c_T \frac{\partial T}{\partial t}\)

• \(\rho\)为材料密度，单位\(kg/m^3\)
• \(c_T\)为材料的比热单位\(J/(kg\cdot K)\)
• \(\kappa _x,\kappa _y,\kappa _z\)分别为沿\(x,y,z\)方向的热传导系数，单位\(W/(m\cdot K)\)
• \(Q(x,y,z,t)\)为物体内部的热源强度，单位\(W/kg\)

3、传热问题的三类边界条件

第一类\(BC(S_1)\)，Dirichlet条件，给定温度分布
\(T(x,y,z,t) =\bar T(t) \quad On \quad S_1\)
\(\bar T(t)\)为在边界\(S_1\)上给定的温度

第二类\(BC(S_2)\)，给定热流密度的Neumann条件
\(\kappa _x \frac{\partial T}{\partial x} n_x+\kappa _y \frac{\partial T}{\partial y} n_y+\kappa _z \frac{\partial T}{\partial z} n_z=\bar q_f \quad On \quad S_2\)
\(n_x,n_y,n_z\)为边界外法线的方向余弦
\(\bar q_f(t)\)为在边界\(S_2\)上的给定热流密度

第三类\(BC(S_3)\)，给定对流换热的Neumann条件
\(\kappa _x \frac{\partial T}{\partial x} n_x+\kappa _y \frac{\partial T}{\partial y} n_y+\kappa _z \frac{\partial T}{\partial z} n_z=\bar h_c(T_\infty -T) \quad On \quad S_3\)
\(h_c\)为在边界\(S_3\)上的对流换热系数\(T_\infty\)为环境温度

物体\(\Omega\)的边界为\(\partial \Omega =S_1+S_2+S_3\)

4、传热问题的变分原理

考虑问题的初始条件IC(initial condition)
\(T(x,y,z,t =0)=\bar T_0(x,y,z)\)

变分提法为在满足三类边界条件及初始条件的许可温度场中，真实的温度场使以下泛函取极小值

\(min _{T\in BC(S_1,S_2,S_3),IC} I=\frac{1}{2} \int _{\Omega}[\kappa _x (\frac{\partial T}{\partial x})^2+\kappa _y (\frac{\partial T}{\partial y} )^2+\kappa _z (\frac{\partial T}{\partial z} )^2- 2(\rho Q-\rho c_T \frac{\partial T}{\partial t})T ]d\Omega\)

在实际问题的处理过程中，第二类和第三类边界条件事先较难满足，因此可将这两个条件耦合进泛函中，即

\(min _{T\in BC(S_1),IC} I=\frac{1}{2} \int _{\Omega}[\kappa _x (\frac{\partial T}{\partial x})^2+\kappa _y (\frac{\partial T}{\partial y} )^2+\kappa _z (\frac{\partial T}{\partial z} )^2- 2(\rho Q-\rho c_T \frac{\partial T}{\partial t})T ]d\Omega\)
\(-\int _{S_2}\bar q_fTdA+\frac{1}{2}\int _{S_3}\bar h_c(T_{\infty} -T)^2dA\)

5、稳态传热问题的有限元分析列式

稳态问题（steady problem）：温度不随时间变化\(\frac{\partial T}{\partial t}=0\)

单元传热矩阵

6、瞬态传热问题的有限元分析列式

温度随时间变化\(\frac{\partial T}{\partial t}\ne 0\)

• 瞬态传热问题求解的特点

• 瞬态传热问题可以转化为一组以时间t为独立变量的线性常微分方程组

影响整个问题的收敛性，必须根据解的稳定性理论来给出一个最大收敛步长的条件，以得到稳定的解，并控制好计算误差

7、平面3节点三角形传热单元

考虑三种情况

1. 无传热边界，即完全为内部单元
2. 如果该单元的\(jm\)边为第二类传热边界\(BC(S_2)\)时：由\(\bar q_f\)参数来描述
3. 如果该单元的\(jm\)边为第三类传热边界\(BC(S_3)\)时：由\(\bar h_c\)参数来描述

7.1、完全为内部单元（无传热边界）

7.2、\(jm\)边为第二类传热边界\(BC(S_2)\)

7.3、\(jm\)边为第三类传热边界\(BC(S_3)\)

传热问题的有限元分析实例

无限长平板稳态温度场的有限元分析

问题描述

模型分析

1、结构的离散化与编号

节点的温度列阵：\(\mathbf{q} _{T}=\begin{bmatrix} T_1 & T_2 & T_3 & T_4 & T_5 & T_6\end{bmatrix} ^T\)

2、单元描述

对于单元3，也是内部单元，与单元2类似，其系节点温度为\(\begin{bmatrix} T_5 & T_4 & T_3\end{bmatrix} ^T\)

3、建立整体有限元分析方程

4、边界条件的处理及方程求解

该问题的边界条件为\(BC(S_3)\)，在前面计算单元的相关矩阵时已作考虑；因此可直接对方程进行求解

理论解

热应力问题的有限元分析：基本原理

1、热应力问题的物理方程

物理方程

指标形式

2、热应力问题的虚功原理

热应力问题中只有物理方程发生了变化

3、热应力问题的有限元分析列式

单元的节点位移列阵

单元的三大类变量场

取虚位移、虚应变的变分

代入虚功原理

热应力问题的有限元分析实例

一个夹持杆结构的温度应力分析

• （1）在\(20^{\circ}C\)时，一个力作用在一个夹持杆件上\(F=300\times 10^3N\)
• （2）随后温度上到\(60^{\circ}C\)，在这种情况下，进行该结构的有限元分析

1、结构的离散化与编号

离散成两个杆单元

2、单元的描述

3、建立整体刚度方程

组装单元，整体刚度矩阵\(\mathbf K \mathbf q=\mathbf P\)

4、边界条件的处理及刚度方程求解

对于位移边界条件，即自由度1、3固定

解得，

5、计算单元的应力

14专题 有限元分析的典型

project Case study A small project with help of the finite element analysis

基本建模Project1：2D问题带孔平板的受力分析

建模要点：
1、属于平面应力问题，单元选择：PLANE182
2、分别建立平面方板和圆孔平面
3、通过布尔运算生成带孔平板
4、对几何模型的线设置网格大小后进行网格划分

APDL：通过日志文件，可以查看apdl命令流记录的操作，具有累赘多余GUI的垃圾操作记录

%%%%%%%%[基本建模Project1]%%%% begin %%%%%%
/PREP7  !pre-processor
ET,1,PLANE182    !selectelementtype(no.1plane182)
KEYOPT,1,3,3   !set plane stress with thickness
R,1,1,    !realconstant(thickness=1)
UIMP,1,EX,,,2.1e5,   !elastic modulus
UIMP,1,PRXY,,,0.3,    !poission ratio
BLC4,0,0,100,100   !create a rectanglar area (x=0,y=0, width=100,height=100), area No.1
CYL4,50,50,5   !create a circular area (center x=50,y=50,rad=5),area No.2
ASBA，1，2    !subtract area No.2 from area No.1, i.e. the area No.1 - the area No.2
ESIZE,0,5    !divide5piecesfor every line
MSHAPE,0,2D    !key=0 for quadrilateral-shaped element(2D)
MSHKEY,0    !free meshing (0)
AMESH,all    !meshall area
FINISH    !pre-processor end
/SOLU    !enter solution environment (for DOF constraints,force,solve)
NSEL,S,LOC,X,0      !select the nodes at x=0
D,all,ux   
NSEL,R,LOC,Y,0     !apply ux=o for selected nodes
D,ALL,UY    !re-select the node at y=o based on above selection (x=o, y=0)
LSEL,S,LOC,X,100      !apply uy=o for selected node Iselect the line at x=0
SFL,all,PRES,-100  !apply a pressure on selected line
ALLSEL    !select all 
SOLVE    !solve
FINISH    !end the solution
/POST1   !entersolutionenvironment(for DOF constraints,force,solve)
PLNSOL,S,X  !displaythedistributionof gxx
!%%%%%%%%[基本建模Project1]%%%%end%%%%%%

（1）如果希望将方板的宽度和高度设为参数（每个变量不超过8个字符）：
plate_w=80
plate_h=120
(2）如果希望将中间孔的位置和半径设为参数：
hole_x=30
hole_y=40
hole_r=8
(3)将弹性模量设为参数
e_modu=1e5
(4)将每边的单元分段设为参数
line_div=6
(5)将外载值设为参数
pressure=200

将参数设置

基本建模Project2：3D问题花形卡盘的扭转受力分析应用建模

花形卡盘的扭转受力分析命令流

Project3：振动模态分析斜拉桥的模态分析

一般的模态分析流程

应用建模Project4：弹塑性分析厚壁圆筒受内压的弹塑性分析应用建模

如图所示厚壁圆筒内壁受到\(120MPa\)均布压力作用的，其内半径\(a=15mm\)，外半径\(b=30mm\)。圆筒材料是理想弹塑性材料，弹性模量\(200GPa\)，泪松比\(0.3\)，屈服极限\(200MPa\)。材料服从von Mises屈服条件。当圆筒的长度很长时，试分析厚壁圆简受力后的弹塑性变形。

建模要点

1. 由于圆筒很长，因此受均布压力的圆简属于平面应变问题，根据对称性用1/4横截面建立计算模型。在ANSYS环境中设置分析类型，选择平面单元PLANE183，设置为平面应变
2. 通过命令< TB,BKIN >将材料设定为双线性动态强化的材料模型，通过命令< TBDATA >给出材料的屈服应力以及切向模量，对理想弹塑性材料，可以将切尚模量设置为较小的数，这里设置为零
3. 对于非线性的弹塑性分析，通过命令< TIME >定义加载步和计算步，采用命令< DELTIM >设定计算子步的步长。分析类型设定为静态分析
4. 在一般后处理中，通过命令< RSYS,1 >设置成柱坐标系，通过命令< PATH >及< PPATH >来定义中间截面的路径，采用命令< PDEF >将应力映射到该路径上，采用命令< PLPATH >显示计算结果
5. 在时间历程后处理中，采用命令< ANSOL >定义平均节点应力作为变量采用命令< PLVAR >画出变量随时间(外载)的变化曲线。

Project5：传热分析钢制圆柱冷却过程的瞬态传热分析应用建模

如图所示，直径为30mm、高度为60mm的钢制实心圆柱。圆柱的初始温度为\(T_0=600\)℃，为均匀分布，对该物体进行热处理，即将圆柱置于温度为\(T_f=20\)℃、换热系数为\(h=1200W/(m^2\)℃）的介质中冷却。

将钢材的密度设为\(7840kg/m^3\)，并忽略密度随温度的变化。这里考虑钢材的热传导系数随温度变化的情况，在0至600℃范围内 \(k=45-0.03T\)（单位\(W/(m\)℃））；比热也随温度变化，在0至600℃范围内，\(C_p=550+0.147T\)（单位\(J/(kg\)℃)）。

计算圆柱冷却过程的温度分布及温度随时间的变化冷却时间为100秒

建模要点】

1. 根据热传导的对称性，此问题可以作为轴对称问题进行分析，取圆柱纵截面的1/4作为计算模型。在圆柱的外侧面和上端面存在对流换热。建立几何模型时长度单位取mm，统一物理量的单位，导热系数的单位取\(W/(mm\)℃)，密度的取\(kg/mm^3\)，换热系数的单位取\(W/(mm^2\)℃)。
2. 采用传热单元plane77，设置材料导热系数和热容随温度线性变化，
3. 设定为瞬态计算< ANTYPE,TRANS >，通过命令< TUNIF >使得初始温度，在圆柱的侧面和端面上定义对流换热边界条件，通过命令< TIME>设定瞬态过程的计算时间
4. 在一般后处理< /POST1>中，通过命令< /EXPAND >进行对称映射显示设置，以云纹图或等值线方式显示温度分布，以失量图方式显示热流分布.在时间后处理< /POST26>中通过命令< NSOL>设置圆柱中心、侧表面中间、端面与侧表面交界处等位置节点上温度变量随时间的变化，在通过命令< PLVAR >显示温度变量的曲线，

Project6：热应力分析桁架结构的温度及装配应力分析

如图所示为一个桁架结构，分析下列两种情形下的节点位移和单元应力：

• a）构件1、3、7和8的温度升高50℃
• b）由于制造误差，构件9和10短了\(0.63mm\)，而构件6长了\(0.27mm\)，但必须进行强制装配架所用材料的相关参数：弹性模量\(E=200GPa\)，线膨胀系数\(α=1.25×10^{-5}(/\)℃)
• c）每个构件的横截面积均为\(100mm^2\)

建模要点

1. 选择可以施加温度载荷的一二维杆单元LINK1，温度可以作为体力施加到单元上。情形（a）是单纯的温度应力问题，情况（b）属于装配应力问题，也可以转化成温度应力问题。
2. 对于情形（a），通过命令< TREF >把参考温度设为20℃，通过命令< BFE >对单元1、3、7和8施加温度70℃，
3. 对于情形（b），同样可以采用命令< BFE >对单元6、9和10上施加适当的温度变化值，使得各自的温度收缩直等价于制造误差，即假定单元由温度变化引起的自由伸循量，与单元长度的差异量相同，由\(\Delta T=\Delta L/(α\cdot t)\) 计算得到等效温度的变化。单元6长了0.27mm，相当于把单元6的温度升高24℃。单元9和10短了0.63mm，相当于把单元9和10的温度降低39.6℃。

高级建模Project7：结构的概率分析——大型液压机机架的概率设计分析

随机输入参数对系统的影响

如图所示，一个大型液压机机架，在其外圈承受有预紧压力，在内圈承受有工作压力，假设由于制造的原因，使得机架的尺寸有一定的误差，同时，所承受的载荷也有一定的分散度，因此，可以将这几个参数考虑为具有概率分布的随机变量，然后，对机架的最大等效应力(MAXVON)、应力强度(MAXINT)、第一主应力(MAX1)进行分析，相关的参数如下

• 材料：\(E=2.1×1011Pa\)，\(\mu=0.3\)
• 几何：\(H=18m\)，\(D=3.4m\)；
  R1服从Gauss分布，均值为2.25m，方差为0.1；
  R2服从Gauss分布，均值为4.5m，方差为0.2；
• 载荷：PRS1服从Gauss分布，均值为261.4MPa，方差为26.07；
  PRS2服从Gauss分布，均值为196.1MPa，方差为19.8；

结构概率分析的一般流程

1. 首先构建宏文件，以便进行概率分析的循环运行，将常规的有限元建模、分析、后处理以及参数提取都放到宏文件中

2. 调入并执行宏命令

3. 进入概率设计分析模块，进行所有与概率设计分析有关的操作包括设置自变量的概率分布、响应参数以及相关系数，并执行概率设计计算，观察结果，图示或列出变量的内容等。

建模要点

• 1.采用单元PLANE182，在几何建模时需要考虑各个几何参数之间的匹配，因为结构的概率分析需要自动改变各个参数以进行多次的参数化儿模型修改和单元战划分；
• 2.采用命令< *CREATE>建立宏文件，给出完整的有限元建模过程，并设置针对随机变量分析的参数，在后处理中，采用命令<*GET>获取需要进行随机响应变量分析的参数；以便于在概率分析时进行循环运行；

3. 在宏文件中需要对将进行概率设计的输入随机变量事先设置为参数，在宏文件的后处理中，还要设置为响应随机变量的参数，采用命令< *GET>进行提取，并设置成为参数；
4. 在主程序中，首先调入宏文件，通过命令< /PDS>进入概率设计分析模块，通过命令< PDANL>指定宏文件进行概率设计，通过命令< PDVAR>设定输入及响应随机变量，对于输入随机变量，指定概率分布的形式及参数;；
5. 通过命令< PDMETH>设置随机处理的方法，由命令< PDLHS>设置模拟循环次数，由命令< PDEXE>运行概率设计分析；
6. 在概率设计分析模块中，直接采用命令< PDHIST>、< PDSHIS>、< PDSENS>、< PDSCAT>给出随机变量的样本分布、采样曲线、概率敏感图、两个相关变量相关性的图形，采用命令< *GET>获取随机变量的平均值及均方差值等信息。

高级建模Project8：p方法的建模与应用——平面问题的p型单元建模与分析

• p方法是保持有限元的剖分网格固定不变，增加各单元上基底函数的阶次，使分析结果逐步向正确解逼近。
• p方法优点：在使用者没有定义精细网格的情况下，仍可获得具有较高精度的解，这些是传统的h方法所做不到的。
• P方法中基底函数的阶次越高，对正确解的逼近程度就越好，但是由于计算机容量和速度的限制，再加上高阶函数会出现插值波动现象，一般情况下多项式函数的最高阶次p<9
• p方法的收敛准则可以包括计算模型节点上的位移应力和应变或整体的应变能

如图所示，带孔方板受均布拉力作用，使用p型单元，分析孔边缘最大应力，相关的参数如下。

• 几何：L=150mm，H=80mm，R1=5mm，R2=10mm，t=0.5mm
• 材料：E=2.1e5Mpa，μ=0.3;
• 载荷：Pressure=100MPa

p方法分析的一般流程

1. 前处理
2. 设置p方法求解，给出p型单元的阶次范围
3. 建立模型，进行网格设置及划分
4. 设置p型单元的误差及收敛参数
5. 在后处理中考察结果

【建模要点】

1. 采用p型单元PLANE145，采用常规的建模流程进行建模，可以采取较粗的网格，p方法将采用增加各单元基底函数阶次的方法来改善计算精度；
2. 根据对称性，取一半的对象进行常规的建模。在该模型中，将生成两个圆面，可以采用命令< WPOFFS>来平移工作面，以便采用命令 < PCIRC>来生成圆面；
3. 在求解之前，需要设置关键节点的物理量，通过命令< PCONV>来采用p方法中的p阶次来以控制计算误差，一般情况下，对局部的计算精度来进行控制，为提高计算效率，对一些关键点周边的单元，可以保持单元的p阶次不变，采用命令< PMOPTS>来进行设置可以设置多点的p方法收敛准则；
4. 求解后，在后处理中，需要采用命令< SET>调出计算结果，采用命令< *GET>获取关键位置上的计算结果，采用命令< PLCONV>及< PPLOT>显示和图示p方法的收敛曲线及单元的p阶次。