计算力学(九)方程组求解和后处理
线性方程组的求解
求解线性方程组的常见方法
- 高斯消元法
- LU分解法
- Cholesky分解法
- 高斯-赛德尔迭代
- QR分解法
- 平方根法
具体细节看这个答案,在最后给出代码
施加载荷
集中力
根据有限元平衡方程\(\mathbf{KV}=\mathbf{P}\\\)直接更改\(P\)即可
体积力
设单位体积内的体积力为 \[ \mathbf f=[f_x,f_y]^T\\ \] 体积力在单元上做功 \[ \displaystyle\iiint_{\Omega_e}\mathbf{u}^T\mathbf{f}\mathrm d\Omega=\mathbf{V}^T_e\displaystyle\iiint_{\Omega_e}\mathbf{N}^T\mathbf{f}\mathrm d\Omega=\mathbf{V}^T_e\mathbf{P}^B_e\\ \begin{align*} \mathbf{P}^B_e&=\displaystyle\iiint_{\Omega_e}\mathbf{N}^T\mathbf{f}\mathrm d\Omega\\ &=\displaystyle\iiint_{\Omega_e}[N_1f_x\quad N_1f_y\quad \cdots\quad N_8f_x\quad N_8f_y]\mathrm d\Omega\\ \end{align*}\\ \] \(f_x,f_y\)就是密度乘以对应方向上的加速度,上式表示将体积力转换为节点等效载荷.
表面力
与体积力同理,表面力也转换为节点等效载荷 \[ \begin{align*} P_e^q&=\displaystyle\iint_{\Gamma_e}\mathbf{N}^T\mathbf{q}\mathrm d\Gamma\\ &=\displaystyle\iint_{\Gamma_e}[N_1q_x\quad N_1q_y\quad \cdots\quad N_8q_x\quad N_8q_y]\mathrm d\Omega\\ \end{align*}\\ \]
单元应力计算
\[ \sigma=[\sigma_x,\sigma_y,\tau_{xy},\sigma_\theta]=\mathbf{D\varepsilon}\\ \]
由于节点所在的每个单元中都会计算节点应力,并且加到对应程序储存的应力处,因此在实际程序中需要除以节点应力被计算的次数,最终得到的是节点平均应力.
模型/约束及载荷的输入
这里给出我用的输入文件约定
- 第一行是模型基本信息,从左到右依次排列,行内所有数据用空格隔开:
- 节点总数
- 单元总数
- 单元维数
- 杨氏模量
- 泊松比
- 单元厚度
- 模型状态:=0对应于平面应力情况,=1对应于平面应变情况,=2对应于轴对称情况
- 高斯积分点数
- 第二行开始时节点编号和节点坐标,每行总共3个数据
- 节点数据结束后是单元编号和单元节点编号节点顺序如下图所示,每行9个数据,支持矩形单元和三角形单元混用,三角形单元使用时最后两个数据为-1
- 体积力,第一行为0时,代表无体积力,为1时代表有体积力,接下来分别输入(1行3个数据)
- 密度
- \(x\)方向加速度
- \(y\)方向加速度
- 集中力,第一行为0时,代表无集中力,为其他正整数对应集中力数量,接下来分别输入(每行3个数据)
- 集中力作用的节点编号
- 该节点的自由度(代表这个集中力分量的方向)
- 集中力分量的大小
- 表面力,第一行为0时,代表无集中力,为其他正整数对应集中力数量,接下来分别输入(每行6个数据)
- 表面力作用在三个节点上的大小
- 三个节点的节点编号
- 约束,第一行对应约束数量,接下来分别输入(每行3个数据)
- 节点编号
- 节点自由度(代表这个约束的方向)
- 位移大小
实际的输入文件在我的github中可找到
结果可视化
利用tecplot进行结果的可视化,tecplot的输入文件有格式要求,具体编码格式如下
VARIABLES="x" "y" "xxxxx" "xxxxx"
这代表接下来输入数据的变量名称,前两个是节点坐标,后面"xxxxx"是节点上数据的名称,可以有多个数据
ZONE NODES=26,ELEMENTS=2.3000000e+01,DATAPACKING=POINT,ZONETYPE=FEQUADRILATERAL
前两个代表节点数量和单元数量,最后一个指用四边形编码方法.
接下来就是按照上面变量名输入数据
最后是输入可视化的单元,对于二维情况tecplot只能输入4个节点的单元,如果单元只有三个节点,就把最后节点重复一次.可能会使用到的6节点三角单元和8节点矩形单元就需要进行分解,这里按个人喜好分解就行.
最终结果
- 主程序按照约定格式读入输入文件中的信息
- 有限元程序求解
- 输出tecplot格式的dat文件
- 在tecplot中可视化
这样最终就可以看到我们有限元程序的结果了!这里计算悬臂梁受集中力载荷的情况,分别用纯矩形单元计算,和矩形单元三角单元混用计算.
可以看到两种情况有点差别,主要由于划分的单元太少导致,程序是没有问题的.
至此,主要内容就完结,撒花!
之后可能会补充一次,三维的整理好了也可能会再补充一次.写这9篇把之前学的又过了一遍,收获不少,之后还是要勤做整理.回头看来,其实都是很简单基础的东西,写的也很简要,但也断断续续写了好久.最近看了看deal.ii,还有很多要学习的.
祝看到这里的朋友们,生活学习科研顺利!感谢!
附上这个专栏的github仓库,不过可能得等两天才能把代码和输入输出文件补充完整.
#! https://zhuanlan.zhihu.com/p/698204031