# How to define RHS in UEL for a static analysis?

Hi every one!

I'm working on developing a shell element to be implemented in ABAQUS through UEL.

It's going through a static analysis. and as you may know for static analysis, ABAQUS just needs AMATRX and RHS to be defined via UEL. It must be mentioned that the stiffness matrix is changing at each time increment. The force that has been applied, is a concentrated force and a function of time.

I've defined the RHS as:

RHS=-[K]*{U}

where [K] is the stiffness (AMATRX) matrix and {u} is the displacement vector at each time increment.

Since I have some difficulties in getting the expected answer, I'm not sure about correctness of this RHS (I've checked B.C., stiffness matrix and etc. and all of them are correct for sure).

I would appreciate if you could help me with defining this RHS.

Best Regards,

### Hi, I suppose you'll get

Hi, I suppose you'll get expected answer for case of using K * U, if you use the constitutive equqtions: stress = Elastic Tensor * B matrix * displacement.

However, if you use RHS =∫[BT]{σ}dv, you'll always get the right answer. Because usually we use return mapping method to calculate the stress, which doesn't equal to the above 'Elastic Tensor * B matrix * displacement', for nonlinear constitutive  equations.

You have mentioned that RHS=∫[BT]{σ}dv the answer is correct for all of the cases. Would you please let me know how I can calculate stress in this case to substitute in the integration?

Regards,

### Hi ashkan,

Hi ashkan,

Abaqus passes displacement and its increments into subroutine UEL.

Then, you may uses kinematic relations to obtain strain from the given displacement and incremental displacement.

Once you gain the strain, the stress would be obtained from the constitutive equations (It depends on the constitutive model you choose).

For linear elastic material, the stress equals to [Elastic tensor] * {Total Strain};

For elastoplastic material, the stress should be updated from the mostly used SIMO's return mapping method (You may find it from the book "[J.C.SIMO]Computational Inelasticity". You can also find the 1D case from the link wrote by yawlou )

### Hi, all I want your help for my UEL subroutine

Actually i have written UEL subroutine  for 4 noded quarilateral element for cantilever beam having end point load. I havewritten input file also.

But it is showing running but not getting run. and gives warning as

***WARNING: 40 elements have incorrect property definitions. The elements have

been identified in element set WarnElemIncorrectProperty.

UEL suroutine:-

SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,

&  PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME,

C

INCLUDE 'ABA_PARAM.INC'

C

DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),PROPS(*),

&  SVARS(*),ENERGY(8),COORDS(MCRD,NNODE),U(NDOFEL),

&  DU(MLVARX,*),V(NDOFEL),A(NDOFEL),TIME(2),PARAMS(*),

&  PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*)

c Variables defined in the UEL subroutine

c   RHS   : Right-Hand-Side vector

c   AMATRX: Stiffness (Jacobian) matrix

c   LFLAGS : Array containing the flags that define the current solution

c            procedures and requirements for element calculations.

c   PREDEF : Array containing the values of predefined variables.

c

c Variables updated in the UEL subroutine

c   SVARS : Maximum separation at each integration point

c

c Variables available in the UEL subroutine

c   U     : Nodal displacement

c   COORDS: Nodal coordinates of an element

c   PROPS : Parameters from an input file

c   MCRD  : Largest active degrees of freedom

c   NNODE : Number of nodes per element

c   MLVARX : Dimensioning parameter used when several displacements or rightside vectors are used.

c   NPREDF : No. of predefined field variables, including temperature

c            For user elements, Abaqus/Standard uses ONE value for each field variable per node.

c   MDLOAD : Total no. of distributed loads defined on this element

c

PARAMETER ( ZERO = 0.D0, HALF = 0.5D0, ONE = 1.D0, TWO = 2.D0,

&             THREE = 1.D0, FOUR= 4.D0, SIX = 6.D0, EIGHT = 8.D0,

&             GAUSSCOORD = 0.577350269D0 )

C

PARAMETER ( NDIM = 2, NDOF_NODE = 2, NGP = 5 )

C

C      NDIM ....      NUMBER OF SPATIAL DIMENSIONS

C      NDOF_NODE .... NUMBEROF DEGREES OF FREEDOM PER NODE

C      NNODE .....    NUMBER OF NODES PER ELEMENT

C      NGP ....       NUMBER OF INTEGRATION POINTS OR GAUSS POINTS

C

DIMENSION B_m(3,8),tB_D(8,3),t_B_D(8,3),D(3,3),xiv(1,5),

&           etav(1,5),weight_gp(1,5),J_m(2,2),

&           SHAPEN(NNODE), dSHAPEN(NDIM,NNODE),GAUSS_POINTS(2,5),

&           INV_J(2,2), ddSHAPEN(NDIM,NNODE),COORDS_IP(2),

&           K_GP(NDOFEL,NDOFEL), F1(NDOFEL,NDOFEL), t_B(8,3)

C

C      COORDS_IP .... COORDINATES AT INTEGRATION POINTS

C      F1        .... MATMUL(t_B_D , B_m)

C

DOUBLE PRECISION E, NU

C     PARAMETER DEFINED IN INPUT FILE:-

!L_X        = PROPS(1)

!L_Y        = PROPS(2)

!THICKNESS  = PROPS(3)

E          = PROPS(1)

NU         = PROPS(2)

!RHO        = PROPS(6)

C     INITIALISE RHS AND LHS:-

DO K1 = 1, NDOFEL

RHS(K1,1) = ZERO

DO K2 = 1,NDOFEL

AMATRX(K1,K2) = ZERO

END DO

END DO

C      GAUSS OR INTEGRATION POINTS ARRAY:-

GAUSS_POINTS = RESHAPE((/-0.9061798459386640, 0.2369268850561891,

&   -0.5384693101056831, 0.4786286704993665, 0, 0.5688888888888889,

&    0.5384693101056831, 0.4786286704993665, 0.9061798459386640,

&    0.2369268850561891/),(/2,5/))

xiv(1,1) = GAUSS_POINTS(1,1)

xiv(1,2) = GAUSS_POINTS(1,2)

xiv(1,3) = GAUSS_POINTS(1,3)

xiv(1,4) = GAUSS_POINTS(1,4)

xiv(1,5) = GAUSS_POINTS(1,5)

etav(1,1) = GAUSS_POINTS(1,1)

etav(1,2) = GAUSS_POINTS(1,2)

etav(1,3) = GAUSS_POINTS(1,3)

etav(1,4) = GAUSS_POINTS(1,4)

etav(1,5) = GAUSS_POINTS(1,5)

weight_gp(1,1) = GAUSS_POINTS(2,1)

weight_gp(1,2) = GAUSS_POINTS(2,2)

weight_gp(1,3) = GAUSS_POINTS(2,3)

weight_gp(1,4) = GAUSS_POINTS(2,4)

weight_gp(1,5) = GAUSS_POINTS(2,5)

C     LOOP OVER GAUSS POINTS OR INTEGRATION POINTS

DO xgp =1,NGP

xi = xiv(1,xgp)

wt_gpxi= weight_gp(1,xgp)

DO egp =1,NGP

eta= etav(1,egp)

wt_gpet= weight_gp(1,egp)

C           SHAPE FUNCTIONS:-

SHAPEN(1) = (ONE - xi)*(ONE - eta)/FOUR

SHAPEN(2) = (ONE + xi)*(ONE - eta)/FOUR

SHAPEN(3) = (ONE + xi)*(ONE + eta)/FOUR

SHAPEN(4) = (ONE - xi)*(ONE + eta)/FOUR

C          DERIVATIVE OF SHAPE FUNCTIONS:-

dSHAPEN(1,1) = -(1- eta)/FOUR  ! wrt xi

dSHAPEN(1,2) = (1-eta)/FOUR    ! wrt xi

dSHAPEN(1,3) =(1+eta)/FOUR     ! wrt xi

dSHAPEN(1,4) = -(1+eta)/FOUR   ! wrt xi

dSHAPEN(2,1) =-(1- xi)/FOUR  ! wrt eta

dSHAPEN(2,2) = -(1+xi)/FOUR  ! wrt eta

dSHAPEN(2,3) =(1+xi)/FOUR    ! wrt eta

dSHAPEN(2,4) =(1-xi)/FOUR    ! wrt eta.

C          COMPUTE COORDINATES AT THE INTEGRATION OR GAUSS POINTS:-

DO K1=1, 2

COORDS_IP(K1) = ZERO

END DO

DO K1=1,NNODE

DO K2=1,NNODE

COORDS_IP(K2)=COORDS_IP(K2)+SHAPEN(K1)*COORDS(K2,K1)

END DO

END DO

C          FIND JACOBIAN AND ITS INVERSE MATRIX:-

DO i = 1, NDIM

DO j = 1, NDIM

J_m(i,j) = ZERO

INV_J(i,j) = ZERO

END DO

END DO

DO INOD= 1, NNODE

DO IDIM = 1, NDIM

DO JDIM = 1, NDIM

J_m(JDIM,IDIM) = J_m(JDIM,IDIM) +

&                           dSHAPEN(JDIM,INOD)*COORDS(IDIM,INOD)

END DO

END DO

END DO

det_J = J_m(1,1)*J_m(2,2) - J_m(1,2)*J_m(2,1)

IF (det_J  .gt. ZERO) THEN

! jacobian is positive - o.k.

INV_J(1,1) = J_m(2,2)/det_J

INV_J(2,2) = J_m(1,1)/det_J

INV_J(1,2) = -J_m(1,2)/det_J

INV_J(2,1) = -J_m(2,1)/det_J

ELSE

! negative or ZERO jacobian

WRITE(7,*)'WARNING: element',jelem,'has negative Jacobian'

END IF

C           COMPUTE DERIVATIVE OF dSHAPEN:-

ddSHAPEN = MATMUL(inv_J,dSHAPEN)

C           COMPUTE B_m MATRIX:-

B_m = ZERO

DO INOD = 1, NNODE

B_m(1,2*INOD-1)=ddSHAPEN(1,INOD)

B_m(2,2*INOD)=ddSHAPEN(2,INOD)

B_m(3,2*INOD-1)=ddSHAPEN(2,INOD)

B_m(3,2*INOD)=ddSHAPEN(1,INOD)

ENDDO

C           COMPUTE ELASTICITY MATRIX D FOR PLANE STRESS CASE:-

D(1,1)= E/(1-NU**2)

D(1,2)= NU*E/(1-NU**2)

D(1,3)= 0

D(2,1)= NU*E/(1-NU**2)

D(2,2)= E/(1-NU**2)

D(2,3)= 0

D(3,1)= 0

D(3,2)= 0

D(3,3)= ((1-NU)/2)*E/(1-NU**2)

C           COMPUTE KGP:-

t_B=TRANSPOSE(B_m)

t_B_D= MATMUL(t_B, D)

F1= MATMUL(t_B_D , B_m)

K_GP = det_J * wt_gpxi * wt_gpet * THICKNESS * F1

C            ASSEMBLY OF RHS AND LHS:-

AMATRX = AMATRX + K_GP

END DO ! END OF EGP

END DO   ! END OF XGP

RETURN

END

************************************************

And input file is:-

Geometry:

length = 1

height = 0.1

width = 0.01

**

*NODE

1,0,0

2,0.05,0

3,0.1,0

4,0.15,0

5,0.2,0

6,0.25,0

7,0.3,0

8,0.35,0

9,0.4,0

10,0.45,0

11,0.5,0

12,0.55,0

13,0.6,0

14,0.65,0

15,0.7,0

16,0.75,0

17,0.8,0

18,0.85,0

19,0.9,0

20,0.95,0

21,1,0

22,0,0.05

23,0.05,0.05

24,0.1,0.05

25,0.15,0.05

26,0.2,0.05

27,0.25,0.05

28,0.3,0.05

29,0.35,0.05

30,0.4,0.05

31,0.45,0.05

32,0.5,0.05

33,0.55,0.05

34,0.6,0.05

35,0.65,0.05

36,0.7,0.05

37,0.75,0.05

38,0.8,0.05

39,0.85,0.05

40,0.9,0.05

41,0.95,0.05

42,1,0.05

43,0,0.1

44,0.05,0.1

45,0.1,0.1

46,0.15,0.1

47,0.2,0.1

48,0.25,0.1

49,0.3,0.1

50,0.35,0.1

51,0.4,0.1

52,0.45,0.1

53,0.5,0.1

54,0.55,0.1

55,0.6,0.1

56,0.65,0.1

57,0.7,0.1

58,0.75,0.1

59,0.8,0.1

60,0.85,0.1

61,0.9,0.1

62,0.95,0.1

63,1,0.1

*USER ELEMENT, TYPE=U1, NODE=4, COORDINATES=2, PROPERTIES=2, VARIABLES=1

1, 2

*ELEMENT,TYPE=U1,ELSET=CANT_BEAM

1,1,2,23,22

2,2,3,24,23

3,3,4,25,24

4,4,5,26,25

5,5,6,27,26

6,6,7,28,27

7,7,8,29,28

8,8,9,30,29

9,9,10,31,30

10,10,11,32,31

11,11,12,33,32

12,12,13,34,33

13,13,14,35,34

14,14,15,36,35

15,15,16,37,36

16,16,17,38,37

17,17,18,39,38

18,18,19,40,39

19,19,20,41,40

20,20,21,42,41

21,22,23,44,43

22,23,24,45,44

23,24,25,46,45

24,25,26,47,46

25,26,27,48,47

26,27,28,49,48

27,28,29,50,49

28,29,30,51,50

29,30,31,52,51

30,31,32,53,52

31,32,33,54,53

32,33,34,55,54

33,34,35,56,55

34,35,36,57,56

35,36,37,58,57

36,37,38,59,58

37,38,39,60,59

38,39,40,61,60

39,40,41,62,61

40,41,42,63,62

**DEFINE THE MATERIAL

*SOLID SECTION, ELSET=CANT_BEAM, MATERIAL=MILD STEEL

0.01

*MATERIAL, NAME=MILD STEEL

*ELASTIC, TYPE=ISOTROPIC

200E6, 0.3

*Nset, nset=End_node

63

*UEL PROPERTY, ELSET=CANT_BEAM

200E6, 0.3

** SOLVING METHOD

*STEP, Perturbation

*STATIC

*BOUNDARY

1, 1, 2, 0.0

22, 1, 2, 0.0

43, 1, 2, 0.0

63, 2, -1

**DEFINE OUTPUT YOU WISH

*OUTPUT,FIELD,FREQUENCY=1

**OUTPUT OF END_POINT NODE:DISPLACEMENT

*NODE OUTPUT,NSET=End_node

u

** u ---- DISPLACEMENT.

*END STEP

PLease tell me corrections , i will be thankful to u all.

### re: UEL

You've specified the element properties twice - delete *solid section, *material and *elastic lines and it should work. I also notice that you incorrectly defined your parameter THREE = 1.d0 in the UEL.

### mailing list

Subscribe to and seek assistance from the
ABAQUS mailing list
https://groups.yahoo.com/group/ABAQUS
Search the archive before posting.

For an intro to subroutines get the file
http://imechanica.org/files/Writing User Subroutines with ABAQUS.pdf

Good luck

Frank 