首页/文章/ 详情

abaqus文档(混凝土坝动力响应分析)

6天前浏览64

概述

  帖子选自abaqus帮助文档,以一座混凝土坝为例讲解如何用abaqus进行结构动力时程计算、应力分析。官方文档这部分的内容很多,这里重点讲解混凝土坝的地震响应计算(我的本行)。
  算例采用附加质量法考虑地震中坝体与库水的动力相互作用,通过用户子程序UEL将附加质量嵌入到abaqus中。坝体调用abaqus提供的CDP模型来模拟材料非线性。在坝体底部输入两个方向的地震加速度时程。
  帖子的所有代码都可以在官方文档上直接下载,文末也会贴上inp文件和UEL子程序代码。下面是详细的内容。

模型信息

  混凝土坝坝底宽度70    ,坝顶宽度14.8    ,坝体高度103    ,水深91.75    ,上游折坡处宽度为19.25    ,坝体截面详细尺寸为
  采用平面应变(CPS4R)单元离散几何模型,网格图为  采用用户自定义附加质量单元模拟坝体与库水动力相互作用,下图中框起来的部分即为附加质量单元,因为是自定义单元,所以abaqus中只显示了“xx”形状。  坝体材料属性信息为

属性数值单位
Young's modulus31027MPa
Poisson's ratio0.15-
Density2643kg/m³
Dilation angle36.31°
Compressive initial yield stress13.0MPa
Compressive ultimate stress24.1MPa
Tensile failure stress2.9MPa

  在坝体底部输入顺河向和竖直向地震加速度时程曲线,顺河向加速度时程曲线为  竖直向加速度时程时程曲线为  调用abaqus提供的CDP模型来模拟坝体混凝土材料的非线性力学行为,其中混凝土的受拉强化曲线为  混凝土的受拉损伤曲线为  详细的数据在文末的inp文件中。

位移、速度和加速度计算结果

  统计了坝顶和坝踵的相对位移、速度和加速度时程曲线,其中,坝顶坝踵顺河向相对位移曲线为  坝顶坝踵顺河向相对速度曲线为  坝顶坝踵顺河向相对加速度曲线为

坝体材料力学响应

  统计了坝体材料在整个计算过程中的应力、应变和损伤数据,其中,坝踵压应力时程曲线为  10    时刻坝体拉损伤云图为

INP文件

  下面给出计算采用的inp文件,其中包括坝体材料的受拉、受压损伤数据,以及用户自定义附加质量单元的数据,顺河向和竖直向地震加速度数据以外部文件的方式导入,数据量太大,这里没有给出,如需要可去官方文档现在或者私信我。

*HEADING
 KOYNA DAM: SEISMIC ANALYSIS 
 Units - n, m, sec
*********************************************
** Step 1: Gravity load                
** Step 2: Hydrostatic pressure load
** Step 3: Earthquake including hydrodynamic 
**    interactions
**
** Requires FORTRAN subroutine addedmass_uel, 
** and amplitude curves koyna_haccel.inp and 
** koyna_vaccel.inp
**
** Execution command:
** abaqus job=koyna_std user=addedmass_uel
*********************************************
*PREPRINT, MODEL=YES
*NODE
   1 ,  0.00,   0.00
   21, 70.00,   0.00
 2401,  0.00,  66.50
 2421, 19.25,  66.50
 3601,  0.00,  91.75
 3621, 16.17,  91.75
 3801,  0.00, 103.00
 3821, 14.80, 103.00
*NGEN, NSET=NBASE
 1, 21
*NGEN, NSET=NMID
 2401, 2421
*NGEN, NSET=NWL
 3601, 3621
*NGEN, NSET=NTOP
 3801, 3821
*NFILL, BIAS=1.05
 NBASE, NMID, 24, 100
*NFILL, BIAS=0.92
 NMID, NWL, 12, 100
*NFILL
 NWL, NTOP, 2, 100
*NSET, NSET=NOUT
 1, 3801
*ELSET, ELSET=EOUT
 2320, 2420
*ELEMENT, TYPE=CPS4R
 1, 1, 2, 102, 101
*ELGEN, ELSET=DAM
 1, 20, 1, 1, 38, 100, 100
*ELSET, ELSET=WDAM, GENERATE
 1, 3501, 100
*USER ELEMENT, NODES=2, TYPE=U1, PROP=3, COORD=2
 1, 2
*ELEMENT, TYPE=U1
 10001, 1, 101
*ELGEN, ELSET=ADDED_MASS
 10001, 36, 100, 1
*UEL PROPERTY, ELSET=ADDED_MASS
 91.75, 91.75, 1000.0
*SOLID SECTION, ELSET=DAM, MATERIAL=CONCRETE
 1.0,
*MATERIAL, NAME=CONCRETE
*ELASTIC
 3.1027E+10, 0.2
*DENSITY
 2643.0
*DAMPING, BETA=0.00323
*CONCRETE DAMAGED PLASTICITY
 36.31
*CONCRETE COMPRESSION HARDENING
 13.0E+6, 0.0
 24.1E+6, 0.001
*CONCRETE TENSION STIFFENING, TYPE=DISPLACEMENT
 2.9E+6         ,0
 1.94393E+6     ,0.000066185
 1.30305E+6     ,0.00012286
 0.873463E+6    ,0.000173427
 0.5855E+6      ,0.00022019
 0.392472E+6    ,0.000264718
 0.263082E+6    ,0.000308088
 0.176349E+6    ,0.00035105
 0.11821E+6     ,0.000394138
 0.0792388E+6   ,0.000437744
 0.0531154E+6   ,0.000482165
*CONCRETE TENSION DAMAGE, TYPE=DISPLACEMENT
 0          ,0
 0.381217   ,0.000066185
 0.617107   ,0.00012286
 0.763072   ,0.000173427
 0.853393   ,0.00022019
 0.909282   ,0.000264718
 0.943865   ,0.000308088
 0.965265   ,0.00035105
 0.978506   ,0.000394138
 0.9867     ,0.000437744
 0.99177    ,0.000482165
*BOUNDARY
 NBASE, 1,2
*AMPLITUDE, NAME=HAMP, INPUT=hamp.inp
*AMPLITUDE, NAME=VAMP, INPUT=hvmp.inp
**
*STEP, NLGEOM, UNSYMM=YES
 STEP 1 - GRAVITY LOAD
*STATIC
 1.0E-10, 1.0E-10
*DLOAD
 DAM, GRAV, 9.81, 0, -1
*OUTPUT, FIELD, VARIABLE=PRESELECT, FREQ=10
*ELEMENT OUTPUT
 S,PE,LE,PEEQ,PEEQT,DAMAGEC,DAMAGET,SDEG
*OUTPUT, HISTORY, VARIABLE=PRESELECT, FREQ=1
*ELEMENT OUTPUT, ELSET=EOUT
 SP1
*NODE OUTPUT, NSET=NOUT
 U
*END STEP
**
*STEP, NLGEOM, UNSYMM=YES
 STEP 2 - HYDROSTATIC LOAD
*STATIC
 1.0E-10, 1.0E-10
*DLOAD
 WDAM, HP4, 900067.6, 91.75, 0
*END STEP
**
*STEP, INC=2000, NLGEOM, UNSYMM=YES
 STEP 3 - EARTHQUAKE
*DYNAMIC, HAFTOL=1.0E7
 0.02, 10.0, 1E-15,0.02
*CONTROLS,PARAMETERS=FIELD
1.E-5
*BOUNDARY,TYPE=ACCELERATION,AMPLITUDE=HAMP
 NBASE,1,1,9.81
*BOUNDARY,TYPE=ACCELERATION,AMPLITUDE=VAMP
 NBASE,2,2,9.81
***CONTROLS, ANALYSIS=DISCONTINUOUS
*END STEP

附加质量单元的UEL的程序

  下面是UEL程序,可以直接拷贝过去用,如果不知道怎么用,可以私信我。

      SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,
     1     NDOFEL,NRHS,NSVARS,PROPS,NPROPS,COORDS,
     2     MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME,
     3     KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,
     4     ADLMAG,PREDEF,NPREDF,LFLAGS,MLVARX,
     5     DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,
     6     PERIOD)
C     
      INCLUDE 'ABA_PARAM.INC'
      PARAMETER ( ZERO = 0.D0, HALF = 0.5D0, 
     1     ONE = 1.D0, SEVEN=7.0D0, EIGHT=8.0D0 )
C
      DIMENSION RHS(MLVARX,*),
     1     AMATRX(NDOFEL,NDOFEL),
     2     SVARS(NSVARS),ENERGY(8),PROPS(*),
     3     COORDS(MCRD,NNODE),U(NDOFEL),
     4     DU(MLVARX,*),V(NDOFEL),A(NDOFEL),
     5     TIME(2),PARAMS(3),JDLTYP(MDLOAD,*),
     6     ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),
     7     PREDEF(2,NPREDF,NNODE),LFLAGS(*),
     8     JPROPS(*)
C
C     USER ELEMENT TO MODEL THE HYDRODYNAMIC 
C     PRESSURE OF RESERVOIR ON DAM USING 
C     WESTERGAARD ADDED MASS TECHNIQUE. 
C     
C     PROPERTIES:
C           PROPS(1) = Depth of reservoir
C           PROPS(2) = Y coordinate of water level
C           PROPS(3) = Density of water
C
      RESERVOIR_DEPTH    = PROPS(1)
      YCOORD_WATER = PROPS(2)
      WATER_DENSITY  = PROPS(3)
      YCOORD_ELCENTROID = 
     *    HALF * (COORDS(2,2) + COORDS(2,1))
      ELCENTROID_DEPTH = 
     *    YCOORD_WATER - YCOORD_ELCENTROID
      RHO = (SEVEN/EIGHT) * WATER_DENSITY *
     *    SQRT(RESERVOIR_DEPTH * ELCENTROID_DEPTH)
C     
      ALEN = ABS(COORDS(2,2)-COORDS(2,1))
      AM   = HALF*RHO*ALEN
C
      DO K1 = 1, NDOFEL                      
        DO KRHS = 1, NRHS
          RHS(K1,KRHS) = ZERO
        END DO
        DO K2 = 1, NDOFEL
          AMATRX(K2,K1) = ZERO
        END DO
      END DO

      do k1 = 1, nsvars
         svars(k1) = zero
      end do
C
      IF (LFLAGS(3).EQ.1) THEN
C       Normal incrementation
        IF (LFLAGS(1).EQ.11
     *        .OR.LFLAGS(1).EQ.12) THEN
C        *DYNAMIC
          BETA  = PARAMS(2)
          DADU = ONE/(BETA*DTIME**2)
C     
          AMATRX(1,1) = AM*DADU
          AMATRX(3,3) = AM*DADU
          RHS(1,1) = RHS(1,1)-AM*A(1)
          RHS(3,1) = RHS(3,1)-AM*A(3)
          ENERGY(1)= HALF*AM*(V(1)*V(1)+V(3)*V(3))
        END IF
      ELSE IF (LFLAGS(3).EQ.4) THEN
C       Mass matrix
         AMATRX(1,1) = AM
         AMATRX(3,3) = AM
      ELSE IF (LFLAGS(3).EQ.5) THEN
C       Half-step residual calculation
         RHS(1,1) = RHS(1,1)-AM*A(1)
         RHS(3,1) = RHS(3,1)-AM*A(3)
      ELSE IF (LFLAGS(3).EQ.6) THEN
C       Initial acceleration calculation
         AMATRX(1,1) = AM
         AMATRX(3,3) = AM
         ENERGY(1) = HALF*AM*(V(1)*V(1)+V(3)*V(3))
      END IF
C
      RETURN
      END         

来源:有限元先生

ACTAbaqus非线性材料
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-11-03
最近编辑:6天前
外太空土豆儿
硕士 我们穷极一生,究竟在追寻什么?
获赞 2粉丝 1文章 47课程 0
点赞
收藏
作者推荐

newmark法求解运动方程(附matlab代码)

问题引入在数值仿真中,常常涉及到静力和动力的问题。静力问题通常忽略物体的惯性力和阻尼效应。也就是说,系统中的力只与物体的静态形变和外部载荷有关,而不考虑动态过程中速度和加速度带来的影响,因此静力问题可以用方程AX=B表示,这种方程用线性代数里面的内容就可以求解,主要是A矩阵求逆的计算速度决定了方程的求解速度,这部分的资料很多,在此不多说,这个帖子重点讲解动力问题。动力问题极为广泛,动力问题(DynamicProblem)是指在力学和物理学中研究物体或结构在外力作用下随时间变化的运动、变形和应力响应。与静力问题不同,动力问题考虑了系统的惯性力和阻尼效应,这使得其分析需要解决随时间变化的运动方程。动力问题广泛应用于涉及冲击、振动、地震、波动等现象的工程和物理系统中。动力问题的特点考虑惯性力和加速度:在动力问题中,物体或结构的运动和变形不仅受外力和内力的影响,还受到惯性力的作用。根据牛顿第二定律,惯性力与物体的加速度成正比,因此在动力分析中必须考虑加速度的影响。时间依赖性:动力问题的解通常是时间的函数,要求在整个时间域内求解物体或结构的运动状态。随着时间的推移,外力、位移、速度和应力等物理量会发生变化。系统响应的类型:自由振动:不受外部力作用的振动,仅依赖于系统的初始状态和自然特性。强迫振动:在外力作用下,系统的响应随外力变化而变化,通常伴随频率和振幅的变化。冲击和瞬态响应:系统受到瞬时力(如冲击或爆炸)时的快速响应,通常表现为瞬态振动。阻尼效应:动力系统中还可能存在阻尼,阻尼力通常与速度成正比,并会导致系统的运动逐渐衰减。动力问题中必须考虑阻尼效应对系统振动和响应的影响。动力问题的基本方程:动力问题的基本方程是牛顿第二定律或达朗贝尔原理的推广,称为运动方程(EquationofMotion)。一个简单的动力问题可以用以下形式的运动方程来描述:这个方程说明了外力F(t)如何引起系统的位移、速度和加速度的时间响应。这是一个二阶微分方程,求解这个方程的方法大致分为解析法和数值解法。解析解当然是最精确的,但只有简单的动力系统或边界条件较为规则的问题,动力问题可以通过解析方法求解。例如:单自由度系统的自由振动:对于简单的弹簧-质量系统,运动方程可以通过特征值问题解析求解,得到系统的固有频率和振动模态。正弦激励下的强迫振动:当系统受到周期性外力(如正弦波)作用时,可以通过解析方法求解系统的稳态响应。数值解法是对于复杂的几何形状、边界条件或载荷,解析解往往难以得到,因此通常采用数值方法来求解动力问题。常见的数值方法包括:有限元法(FEM):将连续体离散为有限个单元,通过求解离散的动力方程来求解整个结构的动态响应。Abaqus等软件广泛使用有限元方法进行复杂动力学分析。有限差分法(FDM):将动力方程的时间和空间导数进行离散化,用数值方法逐步推进时间步长,求解各个时间点的响应。显式和隐式时间积分法:显式方法:如中心差分法,适合处理瞬态问题和高频振动问题。隐式方法:如Newmark方法,适合求解稳定性要求较高的动力问题。这里就重点介绍newmark法求解运动方程,包括固定位移边界条件的施加。newmark理论及程序newmark法是一种常用的时间积分方法,用于求解结构动力学中的运动方程,尤其是在有限元分析中求解动力问题(如振动、冲击、地震响应等)。该方法由NathanM.Newmark于1959年提出,它通过离散时间步长,将动力方程转化为一系列可以逐步求解的代数方程,从而计算出位移、速度和加速度的时间响应。newmark法的优点在于它既适用于求解线性问题,也适用于求解非线性问题,并且它在保持数值稳定性和精度之间提供了一种灵活的平衡。根据参数的选择,Newmark法可以实现不同的时间积分方案,如条件稳定和无条件稳定。newmark公式资料非常的多,基本随便拎起来一本涉及到动力学的书籍都会给出该方法的详细公式,如下图所示上面的资料来自于朱伯芳院士的《有限单元法原理与应用》第四版的第499页。朱伯芳院士在北京于2024年2月4日逝世,每次我想起朱院士序言中的内容,我都会对朱院士充满敬意。朱院士为我国大江大河治理和水电开发作出了卓著贡献,我在此借这个帖子沉痛悼念朱伯芳院士!深切缅怀朱伯芳院士!朱院士的书是我学习有限元过程中推公式用的最多的一本书,当时我是看了一些国内其他书籍,觉得完蛋了,学不会有限元了,感觉山重水复疑无路了,但是看到朱院士的这本书,顿时豁然开朗,柳暗花明。这里贴一下书的封面,有需要的可以私信我。书归正传,有了newmark公式之后,先不慌着写程序,我们通过自己的理论求解了刚度矩阵、质量矩阵和阻尼矩阵之后,还要添加边界条件,才能消除刚体位移,至于添加边界条件又是一个很大的领域,我涉及不多,这里重点讲解newmark,就以最简单的固定位移边界条件为例。添加固定边界的方法有乘大数法、划零置一法和划行划列法等等,我在这里采用的是行划列法,即假如模型的第7个节点是全固定的,就把刚度矩阵等的第3*7-2:3*7(三维问题中)的行和列都直接删除,然后将矩阵传进newmark法中进行求解,求解完成后,再将自由度还原,方便后处理,这里贴一下主要的matlab程序。function[U,V,A]=newmark(deltat,beta,kk,mm,cc,ndim,fbc,kinc,ndofel,nodeu)%结果初始化uu=zeros(ndofel-ndim*length(nodeu),kinc);%didpvv=zeros(ndofel-ndim*length(nodeu),kinc);%velaa=zeros(ndofel-ndim*length(nodeu),kinc);%acc%%基本参数计算alpha=0.25*(0.5+beta)*(0.5+beta);%newamrk计算参数a0=1/(alpha*deltat*deltat);a1=beta/(alpha*deltat);a2=1/(alpha*deltat);a3=1/2/alpha-1;a4=beta/alpha-1;a5=deltat/2*(beta/alpha-2);a6=deltat*(1-beta);a7=beta*deltat;%%迭代计算ek=a0*mm+a1*cc+kk;%有效刚度矩阵ekk=inv(ek);%有效刚度矩阵求逆,只进行一次fori=2:kincfprintf("kinc:%dtime:%f\n",i,deltat*(i-1));ft=fbc(:,i);u=uu(:,i-1);v=vv(:,i-1);a=aa(:,i-1);%计算i+1时刻的有效荷载ft=fbc(:,i);f=ft+mm*(a0.*u+a2.*v+a3.*a)+cc*(a1.*u+a4.*v+a5.*a);%计算i+1时刻的位移uu(:,i)=ekk*f;%计算i+1时刻的速度和加速度aa(:,i)=a0.*(uu(:,i)-u)-a2.*v-a3.*a;vv(:,i)=v+a6.*a+a7.*(aa(:,i));end%%还原自由度U=zeros(ndofel,kinc);V=U;A=U;fori=1:length(nodeu)node=nodeu(i);U(ndim*node-(ndim-1):ndim*node,:)=4e20;%在位移向量中标记固定边界的位置V(ndim*node-(ndim-1):ndim*node,:)=4e20;%在速度向量中标记固定边界的位置A(ndim*node-(ndim-1):ndim*node,:)=4e20;%在加速度向量中标记固定边界的位置end%还原自由度index=0;fori=1:ndofel/ndimifU(ndim*i-(ndim-1),1)==4e20index=index+1;elseU(ndim*i-(ndim-1):ndim*i,:)=uu(ndim*(i-index)-(ndim-1):ndim*(i-index),:);V(ndim*i-(ndim-1):ndim*i,:)=vv(ndim*(i-index)-(ndim-1):ndim*(i-index),:);A(ndim*i-(ndim-1):ndim*i,:)=aa(ndim*(i-index)-(ndim-1):ndim*(i-index),:);endendend注意传进newmark函数的刚度矩阵等都是已经添加过边界条件的。数值算例设计悬臂梁受动荷载算例,悬臂梁一端固定,另一端受动荷载,统计悬臂端下边界的重点数据作对比,示意图为荷载曲线为计算总时长为10s,增量步长为0.01s,采用六面体网格离散,离散后的网格图为采用C3D8单元,共计5120个六面体单元,节点总数为6273。设置两种工况工况一:采用abaqus计算全部的结果工况二:采用自编newmark程序计算统计了悬臂端下边界中点的位移曲线如下图速度曲线如下图加速度曲线如下图自编newmark程序求解运动方程,与abaqus的计算结果相对比,最终结果一致(也有可能是碰巧对准了),但是依然不影响在这过程中对newmark加深了理解。详细的数据和参数不在此罗列了,有兴趣的可以私信我获取全部数据。点击卡片关注我们来源:有限元先生

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈