《有限元分析与应用》

标签：分析 有限元 frac mathbf cdot 应用 bar partial sigma

第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\)

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\)

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

梁单元的坐标变化

分布力的处理

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

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

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

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

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

轴对称单元

分布力的处理

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

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

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

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

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

参数单元的原理

典型空间问题的MATLAB编程

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

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

节点编号与储存带宽

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

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

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

C0 单元与C1单元

单元的拼片试验

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

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

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

第11讲 高阶及复杂单元

1D 高阶单元

2D 高阶单元

3D 高阶单元

基于薄板理论的弯曲板单元

子结构与超级单元

第12讲 有限元分析的应用领域引论(1)

结构振动的有限元分析：基本原理

结构振动的有限元分析实例

弹塑性问题的有限元分析：基本原理

弹塑性问题的有限元分析：非线性问题求解

第13讲 有限元分析的应用领域引论(2)

传热问题的有限元分析：基本原理

传热问题的有限元分析实例

热应力问题的有限元分析：基本原理

热应力问题的有限元分析实例

专题 有限元分析的典型

project Case study A small project with help of the finite element analysis

标签：分析,有限元,frac,mathbf,cdot,应用,bar,partial,sigma
From： https://www.cnblogs.com/invo/p/18371212