有限元进阶编程——逐个单元法

易公子 2022年9月3日浏览：1539

Element-by-element 这一概念最早可以追溯到1983年，有限元先驱之一的Hughes为解决大规模自由度的有限元模型在求解总刚度方程时，计算机出现的存储空间不足，计算速度下降的问题，将迭代技术运用到刚度方程的求解上，创立了一种基于共轭梯度法的 Element-by-element 算法，该算法避免了整体刚度矩阵的组装，只需对逐个单元的求解，如此一来，大大降低了内存的需求，计算速度也得到了显著提升。

数学角度中的算法的详情请可以查阅有关矩阵分析的教材，也可以点击查看

有限元术的往期推文

有限元术，公众号：有限元术【数值算法】共轭梯度法（二）-预处理共轭梯度法

故事的开始

终于要讲这个知识点了，从发现它到现在已经有快半年的时间了，一开始的一头雾水到现在的“略知一二”，中间有放弃过，因为已经超出了自己的接受范围，但是每次想起来总会去尝试看一看，run一run，慢慢的思想上就接受了它，只因对数值计算的一份“无目的性的”热爱。

有限元中的迭代求解

这一部分内容可能会比较枯燥，可以先大致看一下流程，再结合文末提供的代码查看，再次强调一下，本期分享只是我的学习心得，认识可能比较浅，详情可以对照上面链接的参考书的P55。

整体刚度方程从数学的角度上看，可以将之看作一个线性方程组：

$\left[ \mathbf{K}_m \right] \left\{ \mathbf{U} \right\} =\left\{ \mathbf{F} \right\}$

第一步：

$\left\{ \mathbf{P} \right\} _0=\left\{ \mathbf{R} \right\} _0=\left\{ \mathbf{F} \right\} -\left[ \mathbf{K}_{\mathbf{m}} \right] \left\{ \mathbf{U} \right\} _0$

开始迭代，第 $k$ 步：

$\left\{ \mathbf{Q} \right\} =\sum_{i=1}^{\text{单元数}}{\left[ \mathbf{k}_{\mathbf{m}} \right] _i\left\{ \mathbf{p} \right\} _i}$

$\left\{ \mathbf{Q} \right\} _k=\left[ \mathbf{K}_{\mathbf{m}} \right] \left\{ \mathbf{P} \right\} _k$

$\alpha _k=\frac{\left\{ \mathbf{R} \right\} _{k}^{T}\left\{ \mathbf{R} \right\} _k}{\left\{ \mathbf{P} \right\} _{k}^{T}\left\{ \mathbf{Q} \right\} _k}$

$\left\{ \mathbf{U} \right\} _{k+1}=\left\{ \mathbf{U} \right\} _k+\alpha \left\{ \mathbf{Q} \right\} _k$

$\beta _k=\frac{\left\{ \mathbf{R} \right\} _{k+1}^{T}\left\{ \mathbf{R} \right\} _{k+1}}{\left\{ \mathbf{R} \right\} _{k}^{T}\left\{ \mathbf{R} \right\} _k}$

$\left\{ \mathbf{P} \right\} _{k+1}=\left\{ \mathbf{R} \right\} _{k+1}+\beta _k\left\{ \mathbf{P} \right\} _k$

迭代过程中， $k+1$ 步与 $k$ 步的位移差小于设定的精度值，或者迭代次数大于预设的迭代次数上限，迭代停止。

其中， $\left[ \mathbf{k}_{\mathbf{m}} \right] _i$ 是第 $i$ 个单元的单元刚度矩阵， $\left\{ \mathbf{p} \right\} _i$ 是前置矩阵的一部分，按照单元节点顺序进行编号，为方便编程，可以将单元刚度矩阵的对角线元素的倒数作为前置矩阵 $\left\{ \mathbf{P} \right\}$ 。

从以上迭代公式可以看出，求解节点位移过程中完全不需要整体刚度矩阵的组装再求解，通过“逐个单元”的形式逐步求出满足条件的位移。若搭配并行计算，计算速度更加显著提升，感兴趣的小伙伴可以尝试以下，我是还没达到那种程度，望以后可以完全掌握。

相关代码

完整版代码可在公众号：易木木响叮当 后台回复：FEA_Smith，即可自动获取《有限元方法编程》（第5版）课后资源。

资源使用：在Simply Fortran环境运行，打开压缩包后，点击进入code_run\FEA.prj，运行主程序即可，FEA.prj是木木已经测试过的工程文件，可放心使用，调试其他程序时保持工程文件有一个主程序。

!-----invert the preconditioner and get starting loads--
 loads=zero 
 READ(10,*)loaded_nodes,(k,loads(nf(:,k)),i=1,loaded_nodes)
 READ(10,*)fixed_freedoms
 IF(fixed_freedoms/=0)THEN
   ALLOCATE(node(fixed_freedoms),sense(fixed_freedoms),                   &
     value(fixed_freedoms),no(fixed_freedoms),store(fixed_freedoms))
   READ(10,*)(node(i),sense(i),value(i),i=1,fixed_freedoms)
   DO  i=1,fixed_freedoms 
     no(i)=nf(sense(i),node(i)) 
   END DO
   diag_precon(no)=diag_precon(no)+penalty 
   loads(no)=diag_precon(no)*value
   store=diag_precon(no)
 END IF
 CALL mesh_ensi(argv,nlen,g_coord,g_num,element,etype,nf,loads(1:),      &
                nstep,npri,dtim,solid)
 diag_precon(1:)=one/diag_precon(1:) 
 diag_precon(0)=zero 
 d=diag_precon*loads 
 p=d 
 x=zero 
 cg_iters=0
!-------------pcg equation solution--------
 pcg: DO 
   cg_iters=cg_iters+1 
   u=zero
   elements_2: DO iel=1,nels
     g=g_g(:,iel) 
     km=storkm(:,:,iel) 
     u(g)=u(g)+MATMUL(km,p(g)) 
   END DO elements_2
   IF(fixed_freedoms/=0)u(no)=p(no)*store 
   up=DOT_PRODUCT(loads,d)
   alpha=up/DOT_PRODUCT(p,u) 
   xnew=x+p*alpha 
   loads=loads-u*alpha
   d=diag_precon*loads 
   beta=DOT_PRODUCT(loads,d)/up 
   p=d+p*beta
   CALL checon(xnew,x,cg_tol,cg_converged)
   IF(cg_converged.OR.cg_iters==cg_limit)EXIT
 END DO pcg
 WRITE(11,'(A,I5)')" Number of cg iterations to convergence was",cg_iters
 WRITE(11,'(/A)')"  Node   x-disp      y-disp      z-disp" 
 loads=xnew
 DO k=1,nn 
   WRITE(11,'(I6,3E12.4)')k,loads(nf(:,k)) 
 END DO