abaqus文档（混凝土坝动力响应分析）

概述

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

模型信息

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

  在坝体底部输入顺河向和竖直向地震加速度时程曲线，顺河向加速度时程曲线为  竖直向加速度时程时程曲线为  调用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         

来源：有限元先生