abaqus动水附加质量计算子程序——三维用户自定义单元

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,QUAD=0.25D0,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 WATER ON WALL OF TANK BY 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
C---------------WARNING!-------------------
C X--UPSTREAM/DOWNSTEAM DIRECTION
C Y--VERTICAL
C Z--LEFT/RIGHT-BANK DIRECTION
C------------------------------------------
C
RESERVOIR_DEPTH = PROPS(1)
YCOORD_WATER = PROPS(2)
WATER_DENSITY = PROPS(3)
AREA=0.0
RHO=0.0
YCOORD_ELCENTROID=QUAD*
* (coords(2,1)+coords(2,2)
* +coords(2,3)+coords(2,4))
ELCENTROID_DEPTH=YCOORD_WATER-
* YCOORD_ELCENTROID
RHO=(SEVEN/EIGHT) * WATER_DENSITY *
* SQRT(RESERVOIR_DEPTH * ELCENTROID_DEPTH)
C
A122=(coords(1,1)-coords(1,2))**2+
* (coords(2,1)-coords(2,2))**2+
* (coords(3,1)-coords(3,2))**2
A232=(coords(1,2)-coords(1,3))**2+
* (coords(2,2)-coords(2,3))**2+
* (coords(3,2)-coords(3,3))**2
A342=(coords(1,3)-coords(1,4))**2+
* (coords(2,3)-coords(2,4))**2+
* (coords(3,3)-coords(3,4))**2
A412=(coords(1,4)-coords(1,1))**2+
* (coords(2,4)-coords(2,1))**2+
* (coords(3,4)-coords(3,1))**2
A132=(coords(1,1)-coords(1,3))**2+
* (coords(2,1)-coords(2,3))**2+
* (coords(3,1)-coords(3,3))**2
A12=SQRT(A122)
A23=SQRT(A232)
A34=SQRT(A342)
A41=SQRT(A412)
A13=SQRT(A132)
AREA1=QUAD*SQRT((A12+A23+A13)*(A12+A23-A13)
* *(A12+A13-A23)*(A23+A13-A12))
AREA2=QUAD*SQRT((A13+A34+A41)*(A13+A34-A41)
* *(A13+A41-A34)*(A34+A41-A13))
AREA=AREA1+AREA2
AM=RHO*AREA
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
C
DO I=1,15
u(i)=zero
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
RHS(1,1)=RHS(1,1)-AM*A(1)
ENERGY(1)=HALF*AM*(V(1)*V(1))
END IF
ELSE IF(LFLAGS(3).EQ.4) THEN
C Mass matrix
AMATRX(1,1)=AM
ELSE IF(LFLAGS(3).EQ.5) THEN
C Half-step residual calculation
RHS(1,1)=RHS(1,1)-AM*A(1)
ELSE IF(LFLAGS(3).EQ.6) THEN
C Initial acceleration calculation
AMATRX(1,1)=AM
ENERGY(1)=HALF*AM*(V(1)*V(1))
END IF
IF(LFLAGS(1).EQ.41) THEN
C frequency
IF(LFLAGS(3).EQ.2) THEN
C Stiffness matrix
AMATRX(1,1)=0.
END IF
IF(LFLAGS(3).EQ.4) THEN
C Mass matrix
AMATRX(1,1)=AM
END IF
END IF
C
RETURN
END