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,
- ashkan khalili's blog
- Log in or register to post comments
- 17314 reads
Comments
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.
Thanks for your response.
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,
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 "[1997][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.
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
C
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.
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