帖子选自abaqus帮助文档,以一座混凝土坝为例讲解如何用abaqus进行结构动力时程计算、应力分析。官方文档这部分的内容很多,这里重点讲解混凝土坝的地震响应计算(我的本行)。
算例采用附加质量法考虑地震中坝体与库水的动力相互作用,通过用户子程序UEL将附加质量嵌入到abaqus中。坝体调用abaqus提供的CDP模型来模拟材料非线性。在坝体底部输入两个方向的地震加速度时程。
帖子的所有代码都可以在官方文档上直接下载,文末也会贴上inp文件和UEL子程序代码。下面是详细的内容。
混凝土坝坝底宽度70 ,坝顶宽度14.8 ,坝体高度103 ,水深91.75 ,上游折坡处宽度为19.25 ,坝体截面详细尺寸为
采用平面应变(CPS4R)单元离散几何模型,网格图为 采用用户自定义附加质量单元模拟坝体与库水动力相互作用,下图中框起来的部分即为附加质量单元,因为是自定义单元,所以abaqus中只显示了“xx”形状。 坝体材料属性信息为
属性 | 数值 | 单位 |
---|---|---|
Young's modulus | 31027 | MPa |
Poisson's ratio | 0.15 | - |
Density | 2643 | kg/m³ |
Dilation angle | 36.31 | ° |
Compressive initial yield stress | 13.0 | MPa |
Compressive ultimate stress | 24.1 | MPa |
Tensile failure stress | 2.9 | MPa |
在坝体底部输入顺河向和竖直向地震加速度时程曲线,顺河向加速度时程曲线为 竖直向加速度时程时程曲线为 调用abaqus提供的CDP模型来模拟坝体混凝土材料的非线性力学行为,其中混凝土的受拉强化曲线为 混凝土的受拉损伤曲线为 详细的数据在文末的inp文件中。
统计了坝顶和坝踵的相对位移、速度和加速度时程曲线,其中,坝顶坝踵顺河向相对位移曲线为 坝顶坝踵顺河向相对速度曲线为 坝顶坝踵顺河向相对加速度曲线为
统计了坝体材料在整个计算过程中的应力、应变和损伤数据,其中,坝踵压应力时程曲线为 10 时刻坝体拉损伤云图为
下面给出计算采用的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程序,可以直接拷贝过去用,如果不知道怎么用,可以私信我。
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