Element-by-element 这一概念最早可以追溯到1983年,有限元先驱之一的Hughes为解决大规模自由度的有限元模型在求解总刚度方程时,计算机出现的存储空间不足,计算速度下降的问题,将迭代技术运用到刚度方程的求解上,创立了一种基于共轭梯度法的 Element-by-element 算法,该算法避免了整体刚度矩阵的组装,只需对逐个单元的求解,如此一来,大大降低了内存的需求,计算速度也得到了显著提升。
终于要讲这个知识点了,从发现它到现在已经有快半年的时间了,一开始的一头雾水到现在的“略知一二”,中间有放弃过,因为已经超出了自己的接受范围,但是每次想起来总会去尝试看一看,run一run,慢慢的思想上就接受了它,只因对数值计算的一份“无目的性的”热爱。
这一部分内容可能会比较枯燥,可以先大致看一下流程,再结合文末提供的代码查看,再次强调一下,本期分享只是我的学习心得,认识可能比较浅,详情可以对照上面链接的参考书的P55。
整体刚度方程从数学的角度上看,可以将之看作一个线性方程组:
第一步:
开始迭代,第
迭代过程中,
其中,
从以上迭代公式可以看出,求解节点位移过程中完全不需要整体刚度矩阵的组装再求解,通过“逐个单元”的形式逐步求出满足条件的位移。若搭配并行计算,计算速度更加显著提升,感兴趣的小伙伴可以尝试以下,我是还没达到那种程度,望以后可以完全掌握。
完整版代码可在后台回复: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))
ENDDO
diag_precon(no)=diag_precon(no)+penalty
loads(no)=diag_precon(no)*value
store=diag_precon(no)
ENDIF
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))
ENDDO 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
ENDDO 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))
ENDDO
上述程序考虑了已知位移边界条件,使用的“乘大数法”处理。
最后说明:学识有限,对于上文中提及的逐个单元法有限元进阶技术,只能分享出一点点的相关,随着对有限元认识的加深,相信会慢慢掌握它的,到时再做一次分享,到时会更加清楚的给大家讲明白。