首页/文章/ 详情

学习有限元最好的资料:开源程序

1天前浏览10

写在前面

  授人以鱼,不如授人以渔。这是一期关于如何更好的学习有限元的帖子,看完保证会给你学习有限元理论打开一条新思路。
  经常翻有限元书籍的家人们应该都有一种感觉,有限元领域的书籍大致可分为两种:一种是工程领域的专家编写的,另一种是数学领域的专家编写的第一种书籍往往从杆单元一路讲到非线性,另一种是从无穷维空间在有限维空间的投影一路讲到偏微分方程求解, 首先毋庸置疑,这两种书籍都是对有限元理论的诠释,只是站的角度不一样。
  假如你以超越常人的毅力按照两者之一学习完毕,准备大显身手,开始编程工作,你一定会突然发现:无从下手!怎么办?
  这个帖子就是为了解决这个问题而写的,我主要讲两个方面的事情,一个是上面的两个路径如何选择,另一个是如何解决写程序无从下手的问题

新手不要陷入数学推导

  不要按照一件事做完之后再做另一件事的线性思路来学习有限元理论。我在前面讲了两种可能的有限元书籍,无论是哪一种,新手经常会拿着书就开始啃,陷入到数学公式的推导中,我认为新手一开始就硬肝理论是一种非常不明智的行为,这种线性的学习思路会带来极大的挫败感,我的建议是:理论方面的东西,第一遍看的时候浅尝辄止,然后就直接上软件,在实践中学习理论,就像是教员交给我们的一句话:在战争中学习战争
  简单过完一遍理论就上手软件,会给人带来极大的成就感,但是这个成就感的时间是非常的短暂的,随着模型逐渐复杂,软件会频繁报错,这一些错误会让我们手足无措,里面的名词甚至我们都看不懂,更谈不上去调试我们的模型,这时候我们的学习就会进入下一个环节:重点精读理论并在操作软件中深刻理解理论
  假如,软件给我们一个错误:刚度矩阵奇异,方程无法求解,如果不知道理论,这个错误就会显得莫名其妙,但是如果我们学习过有限元理论中关于边界条件部分的内容,我们就会马上反应过来:我是不是忘了加边界条件了,没有加边界条件,总体刚度矩阵求逆的时候行列式直接为零了,这不就是矩阵奇异了么
  再比如,我们第一次看到软件的界面的时候,肯定会充满疑惑:为什么这么设计?为什么要画网格?为什么要给网格赋予单元类型?如果我们带着这样的疑问再去学习有限元理论,目的就会更加的明确!学完之后,你就会发现软件的设计是如此的精妙!
  这阶段的学习,是实践与理论相互交错的,软件给我们报错,我们回去学习理论,然后根据理论再来调试我们的模型,这种实践与理论循环的方式是符合事物的发展规律的

手搓有限元

  当我们对软件的报错不慌不忙,很淡定的就能根据软件的报错就判断出来对应的理论,然后根据理论去软件的相应位置去修改,这时候我们对软件的学习就应该告一段落了。我们要进入下一个更高的阶段:手编程序
  手编程序的开始工作非常重要,我的建议是:千万不要自己一上去就开始写程序。虽然我就是直接上去就写程序了,这也导致我现在有非常多的坏毛病。
  先找到一个开源的、质量比较高的FEM库
  手搓有限元的思路,不是按照我们对理论的理解就直接上去一顿乱造,一顿for循环,无论什么数据都是直接double数组,甚至上高维数组。这种程序写出来了可读性很差、不可扩展,俗称:屎山代码!(我现在就处于这个阶段),一定要找一个开源的FEM库,先读别人的代码架构是怎么设计的、学习别人是怎么实现积分的、学习别人的矩阵是怎么储存的、学习别人的单元是什么数据结构,肯定不是不管三七二十一都按照double数组,如果涉及到了C++语言,尤其要学习别人的指针是怎么用的,这是C++语言的核心

OpenSees介绍

  这里以开源程序OpenSees为例,大致讲解OpenSees里的内容,如果有人需要这个代码,可以私信我
  OpenSees(Open System for Earthquake Engineering Simulation)是一个专门用于地震和结构工程模拟的开源软件框架。以下是关于OpenSees的一些简要介绍:

  1. 开发背景:OpenSees由美国加州大学伯克利分校的太平洋地震工程研究中心(PEER)开发。

  2. 应用领域:OpenSees在地震工程中应用广泛,不仅能进行结构的地震响应分析,还能用于结构的优化设计、损伤评估以及地震风险分析。它主要用于结构和岩土方面的地震反应模拟,可以实现包括静力线弹性分析、静力非线性分析、截面分析、模态分析、pushover拟动力分析、动力线弹性分析和复杂的动力非线性分析等。

  3. 功能特点

    • 开源性:OpenSees是完全开源的,允许用户自由地修改和扩展软件功能。
    • 模块化设计:软件采用模块化设计,便于用户根据需要添加或修改分析模块。
    • 强大的分析能力:支持多种分析类型,包括线性、非线性、静力、动力分析等。
    • 丰富的材料模型:提供了多种材料模型,如混凝土、钢材、土等,以满足不同结构的仿真需求。
    • 用户友好:通过Python或Tcl脚本语言,用户可以轻松地创建和运行模型。
  4. 系统架构与核心功能:OpenSees采用了模块化的设计理念,其核心组件包括分析内核、对象模型以及命令解释器等。这种设计使得OpenSees既能够处理复杂的非线性动力学问题,又具备良好的用户友好性。

  5. 安装与配置:OpenSees支持跨平台运行,无论是Windows、Linux还是Mac OS用户都能够轻松上手。

  6. 材料模型应用:OpenSees提供丰富的材料模型,可以实现对结构的多种非线性响应分析。

  7. 实际工程案例应用:OpenSees在实际工程项目中扮演着重要角色,例如在桥梁抗震加固项目和城市轨道交通建设中,OpenSees被广泛应用于地铁隧道的稳定性评估。

  8. 未来发展与趋势:随着人工智能和云计算技术的深度融合,OpenSees有望在提高建筑物抗灾能力方面发挥更大的作用。

  OpenSees因其灵活性、可扩展性以及开源特性,在地震工程领域展现出巨大的潜力,成为了研究人员和工程师们不可或缺的工具。
  这里面主要是一些C++的头文件和源程序文件,对应于abaqus的C3D8单元,还有BBAR修正的程序,这些程序都是可以直接读的,非常有利我们学习编程!这里我把形函数的程序贴在下面,供大家学习


// $Revision: 1.1 $
// $Date: 2001-07-11 21:54:41 $
// $Source: /usr/local/cvs/OpenSees/SRC/element/brick/shp3d.cpp,v $

// Ed "C++" Love
//
// 3-d isoparametric 8-node element shape function
//

/* 
      Purpose: Compute 3-d isoparametric 8-node element shape
               functions and their derivatives w/r x,y,z

      Inputs:
         ss[3]     - Natural coordinates of point
         xl[3][8]  - Nodal coordinates for element

      Outputs:
         xsj        - Jacobian determinant at point
         shp[4][8]  - Shape functions and derivatives at point
                     shp[0][i] = dN_i/dx
                     shp[1][i] = dN_i/dy
                     shp[2][i] = dN_i/dzc
                     shp[3][i] =  N_i
*/

void  shp3d( const double ss[3],
      double &xsj,
      double shp[4][8],
      const double xl[3][8]   )
{
    int   i, j, k ;

    double rxsj, ap1, am1, ap2, am2, ap3, am3, c1,c2,c3 ;

    static double xs[3][3] ; 
    static double ad[3][3] ;


      //Compute shape functions and their natural coord. derivatives

      ap1 = 1.0 + ss[0] ;
      am1 = 1.0 - ss[0] ;
      ap2 = 1.0 + ss[1] ;
      am2 = 1.0 - ss[1] ;
      ap3 = 1.0 + ss[2] ;
      am3 = 1.0 - ss[2] ;

      //Compute for ( - , - ) values

      c1      = 0.125*am1*am2 ;
      c2      = 0.125*am2*am3 ;
      c3      = 0.125*am1*am3 ;
      shp[0][0] = -c2 ;
      shp[0][1] =  c2 ;
      shp[1][0] = -c3 ;
      shp[1][3] =  c3 ;
      shp[2][0] = -c1 ;
      shp[2][4] =  c1 ;
      shp[3][0] =  c1*am3 ;
      shp[3][4] =  c1*ap3 ;

      //Compute for ( + , + ) values

      c1      = 0.125*ap1*ap2 ;
      c2      = 0.125*ap2*ap3 ;
      c3      = 0.125*ap1*ap3 ;
      shp[0][7] = -c2 ;
      shp[0][6] =  c2 ;
      shp[1][5] = -c3 ;
      shp[1][6] =  c3 ;
      shp[2][2] = -c1 ;
      shp[2][6] =  c1 ;
      shp[3][2] =  c1*am3 ;
      shp[3][6] =  c1*ap3 ;

      //Compute for ( - , + ) values

      c1      = 0.125*am1*ap2 ;
      c2      = 0.125*am2*ap3 ;
      c3      = 0.125*am1*ap3 ;
      shp[0][4] = -c2 ;
      shp[0][5] =  c2 ; 
      shp[1][4] = -c3 ;
      shp[1][7] =  c3 ;
      shp[2][3] = -c1 ;
      shp[2][7] =  c1 ;
      shp[3][3] =  c1*am3 ;
      shp[3][7] =  c1*ap3 ;

      //Compute for ( + , - ) values

      c1      = 0.125*ap1*am2 ;
      c2      = 0.125*ap2*am3 ;
      c3      = 0.125*ap1*am3 ;
      shp[0][3] = -c2 ;
      shp[0][2] =  c2 ;
      shp[1][1] = -c3 ;
      shp[1][2] =  c3 ;
      shp[2][1] = -c1 ;
      shp[2][5] =  c1 ;
      shp[3][1] =  c1*am3 ;
      shp[3][5] =  c1*ap3 ;

      //Compute jacobian transformation

      for ( j=0; j<3; j++ ) {

        xs[j][0] = ( xl[j][1] - xl[j][0] )*shp[0][1]
                 + ( xl[j][2] - xl[j][3] )*shp[0][2]
                 + ( xl[j][5] - xl[j][4] )*shp[0][5]
          + ( xl[j][6] - xl[j][7] )*shp[0][6] ;

        xs[j][1] = ( xl[j][2] - xl[j][1] )*shp[1][2]
                 + ( xl[j][3] - xl[j][0] )*shp[1][3]
                 + ( xl[j][6] - xl[j][5] )*shp[1][6]
          + ( xl[j][7] - xl[j][4] )*shp[1][7] ;

        xs[j][2] = ( xl[j][4] - xl[j][0] )*shp[2][4]
                 + ( xl[j][5] - xl[j][1] )*shp[2][5]
                 + ( xl[j][6] - xl[j][2] )*shp[2][6]
          + ( xl[j][7] - xl[j][3] )*shp[2][7] ;

      } //end for j     


      //Compute adjoint to jacobian

      ad[0][0] = xs[1][1]*xs[2][2] - xs[1][2]*xs[2][1] ;
      ad[0][1] = xs[2][1]*xs[0][2] - xs[2][2]*xs[0][1] ;
      ad[0][2] = xs[0][1]*xs[1][2] - xs[0][2]*xs[1][1] ;

      ad[1][0] = xs[1][2]*xs[2][0] - xs[1][0]*xs[2][2] ;
      ad[1][1] = xs[2][2]*xs[0][0] - xs[2][0]*xs[0][2] ;
      ad[1][2] = xs[0][2]*xs[1][0] - xs[0][0]*xs[1][2] ;

      ad[2][0] = xs[1][0]*xs[2][1] - xs[1][1]*xs[2][0] ;
      ad[2][1] = xs[2][0]*xs[0][1] - xs[2][1]*xs[0][0] ; 
      ad[2][2] = xs[0][0]*xs[1][1] - xs[0][1]*xs[1][0] ;

      //Compute determinant of jacobian

      xsj  = xs[0][0]*ad[0][0] + xs[0][1]*ad[1][0] + xs[0][2]*ad[2][0] ;
      rxsj = 1.0/xsj ;

      //Compute jacobian inverse

      for ( j=0; j<3; j++ ) {

        for ( i=0; i<3; i++ ) 
      xs[i][j] = ad[i][j]*rxsj ;

      } //end for j


      //Compute derivatives with respect to global coords.

      for ( k=0; k<8; k++) {

 c1 = shp[0][k]*xs[0][0] + shp[1][k]*xs[1][0] + shp[2][k]*xs[2][0] ;
        c2 = shp[0][k]*xs[0][1] + shp[1][k]*xs[1][1] + shp[2][k]*xs[2][1] ;
        c3 = shp[0][k]*xs[0][2] + shp[1][k]*xs[1][2] + shp[2][k]*xs[2][2] ;

        shp[0][k] = c1 ;
        shp[1][k] = c2 ;
        shp[2][k] = c3 ;

      } //end for k

   return ;
}

  写的非常的美观!可以为我们所用,但是一定要需要遵循开源协议。
  这部分的学习,永无止境,可以一直学习到自己毕业,或者学习到离开有限元领域。

FEM的源头:偏微分方程数值解

  这就要进入下一部分学习了:数学角度的有限元
  这部分就开始涉及到有限元的源头:数学了,偏微分方程是这部分学习的主旋律,一切都是围绕偏微分方程的求解,你会涉及到:偏微分方程的强形式和弱形式,偏微分方程的弱解和数值解等等,还有加权余量法、泛函分析与变分原理等等,我也只知道这么多了。
  有限元是一门博大精深的领域,底层是数学中的偏微分方程求解,上面是工程中的力学原理,再上层是工程问题的抽象处理,有关于有限元的每一个名词,单拎出来都可以作为一个领域去研究,但凡你在一个小分支有了突破,都足够你威震武林,至于所有的分支都涉及到,那是不可能的(手动狗头)。

总结

  学习有限元不可莽撞,要有条有理。最好是先过一遍理论,第一遍浅尝辄止即可,不要贪多;然后就上手有限元软件,在实践中学习理论,在战争中学习战争,这一阶段要实践与理论并重,理论指导实践;在手动编程的时候,不建议直接开始写,最好是找到一个开源的FEM库,学习程序的架构设计、数据组织和储存方式等等。大部分人的有限元学习之路,到这里就结束了。
  若有余力,研习偏微分方程数值解

来源:有限元先生
SystemAbaqus非线性轨道交通建筑python岩土云计算理论Opensees材料人工智能
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2025-01-06
最近编辑:1天前
外太空土豆儿
博士 我们穷极一生,究竟在追寻什么?
获赞 25粉丝 5文章 68课程 0
点赞
收藏
未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习计划 福利任务
下载APP
联系我们
帮助与反馈