有限元初学

最新推荐文章于 2024-10-12 09:16:26 发布

qq102385

最新推荐文章于 2024-10-12 09:16:26 发布

阅读量401

点赞数

文章标签： matlab 算法笔记课程设计

本文链接：https://blog.csdn.net/qq102385/article/details/130063433

版权

本文详细介绍了如何使用MATLAB进行有限元分析，解决薄板圆孔的应力集中问题。作者首先建立了物理模型和有限元模型，包括网格化、结点编号、刚度矩阵计算、边界条件和外部载荷。然后，讨论了位移边界条件和应力计算方法。最后，通过编程求解并进行了误差分析。

摘要由CSDN通过智能技术生成

写在前面

博主为有限元学习相关技术学习小白，以下内容是博主近期课程思考的总结，如有错误还请海涵QAQ.

有限元模型建立

物理模型建立

首先约定以下假设：
假设1：结构为矩形薄板，则问题划为平面应力问题；
假设2：重力沿Z轴方向，认为平板XY方向无重力分量，题目不考虑重力体力；
假设3：为满足题目“离开边界较远处”要求，约定矩阵最短边长度的一半大于10倍的圆半径，即 $\frac{{{l_{\min }}}}{2} \ge 10a$

假设4：为满足薄板要求，板厚满足 $\frac{t}{{{l_{\min }}}} < 0.2$

分析问题，由于薄板受力与结构的对称性，取板的1/4建立有限元模型。其位移边界条件为： $x = 0$ 边界上， $u = 0$
边界上 $y = 0$ ， $v = 0$

基于以上分析，对问题相关参数赋值，用于后续计算。1/4矩形薄板长边 $l o n 1 = 12 c m$ ，短边 $l o n 2 = 11 c m$ ，长边沿x轴，短边沿y轴。圆孔半径 $r = 1 c m$ ，板厚 $t = 0.2 c m$ 。X向均布载荷集度 ${q_1} = 20Mpa$ ，y向均布载荷集度 ${q_2} = 10Mpa$ 。材料常数弹性模量 $\times {10^5}Mpa$ ，泊松比 $v = 0.3$ 。图1为矩形薄板部分数据。
等效矩形薄板

网格化及结点总体编号

选择三节点三角形单元作为划分矩形薄板的基本单元，为简化圆孔附近的单元划分难度，需对圆孔边界进行一定简化。参考圆弧插值的基本思路，将四分之一圆弧简化为一系列短直线的组合，插值思路如下：
从圆弧起始点 ${x_0},{y_0})$ 向圆弧终止点 ${x_1},{y_1})$ 以步长 $(d x, d y)$ 进行逆时针插值，当插值点 ${x_i},{y_i})$ 在圆内，则 ${y_{i + 1}} = {y_i} + dy$ 否则 ${x_{i + 1}} = {x_i} - dx$ 直到 ${x_i} = {x_1},{y_i} = {y_1}$ 插值停止。

经过插值，总能在有限步内将圆弧转化为一系列短直线的集合。将圆弧转化为一系列段直线后，问题图形可使用一系列排列整齐的小矩形组成，极大的简化了网格划分的难度。
总体上将图形划分为四个部分，区域1 $\in [0,4a],y \in (0,4a)$ 区域2 $\in (4a,lon1],y \in (0,4a)$ 区域3 $\in [0,4a],y \in [4a,lon2]$ 区域4 $\in (4a,lon1],y \in [4a,lon2]$ 各区域内小矩形边长不同，其边长可根据精度要求进行调整。
各小矩形在整体中可使用 i 行， j 列来表示其位置，则小矩形位置唯一。将小矩形沿对角线分为两个三节点三角形单元，记为上三角形a，下三角形b。由于小矩形位置唯一由(i,j) 表示，则三角形单元各节点坐标可使用(i,j) 计算，且对应关系唯一，记对应关系为
$\left[ {\begin{array}{ccccccccccccccc}{{x_i}}&{{y_i}}\\{{x_j}}&{y{}_j}\\{{x_m}}&{{y_m}}\end{array}} \right] = f(i,j)$
将各结点顺序编号（就是贪吃蛇或者从左到右），则各结点总体编号与 (i,j) 对应。
$\left[ {\begin{array}{ccccccccccccccc}{nu{m_i}}\\{nu{m_j}}\\{nu{m_m}}\end{array}} \right] = g(i,j)$
利用上述对应关系，能够建立节点坐标与结点总体编号的对应关系。在单元刚度矩阵的相关计算中，直接总体刚度矩阵进行赋值操作，省去存储节点位移和结点总体编号的过程，同时也免去总体刚度矩阵赋值中不断搜索对应的过程，简化计算。（当然如果你喜欢，建议直接暴力存矩阵，遍历对应即可）

刚度矩阵计算

根据网格化结果，矩形薄板均由三节点三角形单元组成，选择线性位移模式，三角形单元位移插值函数为
$\begin{array}{l}u = {N_i}{u_i} + {N_j}{u_j} + {N_m}{u_m}\\v = {N_i}{v_i} + {N_j}{v_j} + {N_m}{v_m}\end{array}$
其中
${N_i} = \frac{1}{{2A}}({a_i} + {b_i}x + {c_i}y)\begin{array}{ccccccccccccccc}{}&{(i,j,m)}\end{array}$
三角形单元位移为
$\left[ {\begin{array}{ccccccccccccccc}{{N_i}}&{{N_j}}&{{N_m}}\end{array}} \right]{a^e} = N{a^e}$
三角形单元应变为
$\varepsilon = L\left[ {\begin{array}{ccccccccccccccc}{{N_i}}&{{N_j}}&{{N_m}}\end{array}} \right]{a^e} = \left[ {\begin{array}{ccccccccccccccc}{{B_i}}&{{B_j}}&{{B_m}}\end{array}} \right]{a^e} = B{a^e}$
其中
$\left[ {\begin{array}{ccccccccccccccc}{{B_i}}&{{B_j}}&{{B_m}}\end{array}} \right] = \frac{1}{{2A}}\left[ {\begin{array}{ccccccccccccccc}{{b_i}}&0&{{b_j}}&0&{{b_m}}&0\\0&{{c_i}}&0&{{c_j}}&0&{{c_m}}\\{{c_i}}&{{b_i}}&{{c_j}}&{{b_j}}&{{c_m}}&{{b_m}}\end{array}} \right]$
单元应力为
$\sigma = \left( {\begin{array}{ccccccccccccccc}{{\sigma _x}}\\{{\sigma _y}}\\{{\tau _{xy}}}\end{array}} \right) = D\varepsilon = DB{a^e} = S{a^e}$
其中弹性矩阵为
$\frac{{{E_0}}}{{1 - \nu _0^2}}\left[ {\begin{array}{ccccccccccccccc}1&{{v_0}}&0\\{{v_0}}&1&0\\0&0&{\frac{{1 - {v_0}}}{2}}\end{array}} \right]$
对于平面应力问题
${E_0} = E,{\nu _0} = \nu$
最后，三节点三角形单元的单元刚度矩阵为
${K^e} = \int_{{\Omega ^e}} {{B^T}DBtdxdy = } {B^T}DBtA$
其中，A为三节点三角形单元面积。
利用结点总体编号规则，将单元刚度矩阵的元素“对号入座”，叠加到结构矩阵中，得到总体刚度矩阵 K。

外部载荷计算

根据题目假设，在不考虑重力作用下，矩形薄板受到的外力仅有分别作用在x向边界的x向均布载荷集度和作用在y边界y向均布载荷集度。则单元等效结点载荷为
${P^e} = P_S^e = \int_{S_\sigma ^e} {{N^T}TtdS}$
边界处单元受力情况一致，取其中一个单元进行分析，设x向载荷作用于边，以边建立局部坐标系，沿边插值函数可写作
${N_i} = 1 - \frac{s}{l}\begin{array}{ccccccccccccccc}{}&{}\end{array}{N_j} = \frac{s}{l}\begin{array}{ccccccccccccccc}{}&{}\end{array}{N_m} = 0$
边界上的面积力为
${\left[ {\begin{array}{ccccccccccccccc}{{q_1}}&0\end{array}} \right]^T}$
单元等效结点载荷为
${P^e} = \frac{1}{2}{q_1}lt{[\begin{array}{ccccccccccccccc}1&0&1&0&0&0\end{array}]^T}$
同理，取边界处单元，可得其等效结点载荷为
${P^e} = \frac{1}{2}{q_2}lt{[\begin{array}{ccccccccccccccc}0&1&0&1&0&0\end{array}]^T}$
由于网格划为四个部分，上式中需根据区域变化。
将边界各结点受力，按总体编号叠加，得到总体外力向量 P。

位移边界条件

在求解位移场问题时，需要提出足以约束刚体位移的几何边界条件，以消除结构刚度矩阵的奇异性。根据模型简化假设，其位移边界要求为： $x = 0, u = 0$ $y = 0, v = 0$ 求解方程时，采用直接带入法修正方程。将已知的结点位移所对应的总体刚度矩阵对应的行列消去，同时将其对应的外力向量对应的列也消去，重新组合方程。消去过程保留原有结点的顺序。方程求解完成后，根据消去规则，将代求结点位移重新对号入座，再结合已知结点位移，得到系统总体位移向量a 。

应力计算

在求得结点位移后，根据公式求得的单元应变和应力，由于导数运算的结果，在单元交界面上应力不连续，力的边界条件也不能精确满足，因而需要对它作必要的处理。此次计算选用单元平均的方法处理以上问题。在三节点三角形单元中，应力解在单元内是常数。由于应力近似解总是在精确解上下振荡，可取相邻单元应力的平均值作为此两个单元合成的较大四边形单元形心处的应力。精确面积加权平均公式如下。

参考网格划分的特点，选取同一小矩形下的上三角形a和下三角形b作为应力平均四边形，简化编程难度。