User login

Navigation

You are here

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,

Comments

Weijie Liu's picture

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.

Thanks for your response.

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,

Weijie Liu's picture

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 "[1997][J.C.SIMO]Computational Inelasticity". You can also find the 1D case from the link wrote by yawlou )

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.

Please help me to get out of this.

UEL suroutine:-

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

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

     &  KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF,

     &  LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD)

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(*),

     &  JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),

     &  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:-

*Heading

 

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

** LOADS AND BOUNDARY CONDITIONS

*BOUNDARY

 1, 1, 2, 0.0

 22, 1, 2, 0.0

 43, 1, 2, 0.0

*CLOAD

   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.

Thanks in advance.

            

           

Jason Mayeur's picture

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.

Frank Richter's picture

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

Subscribe to Comments for "How to define RHS in UEL for a static analysis?"

More comments

Syndicate

Subscribe to Syndicate