首页/文章/ 详情

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

1年前浏览329

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

上述程序考虑了已知位移边界条件,使用的“乘大数法”处理。

最后说明:学识有限,对于上文中提及的逐个单元法有限元进阶技术,只能分享出一点点的相关,随着对有限元认识的加深,相信会慢慢掌握它的,到时再做一次分享,到时会更加清楚的给大家讲明白。

来源:易木木响叮当
理论通用
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2023-06-01
最近编辑:1年前
易木木响叮当
硕士 有限元爱好者
获赞 224粉丝 284文章 355课程 2
点赞
收藏
未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习计划 福利任务
下载APP
联系我们
帮助与反馈