As per this documentation page UEL - SIMULIA User Assistance 2022, the value of AMATRX is to be set based on the value of LFLAGS variable. In my case, considering matrix generation, it seems that with LFLAGS(3)=2, the stiffness matrix is correctly assembled, but, instead, with LFLAGS(3)=4, the mass matrix is completely null, which shouldn’t be. However, printing the values of the latter before assignment yields the correct results. How should I proceed? In case, I also attach the UEL source code, the subroutine in question is the UPE2, which is associated with planar elements.
The commands that I have used to generate the matrices are
*MATRIX GENERATE, STIFFNESS, MASS=CONSISTENT, SOURCE=ELEMENTS
*MATRIX OUTPUT, STIFFNESS, MASS, FORMAT=MATRIX INPUT
The code in question is
SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,
+ PROPS,NPROPS,coords,MCRD,NNODE,Uall,DUall,Vel,Accn,JTYPE,
+ TIME,DTIME,KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,
+ PREDEF,NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,
+ NJPROP,PERIOD)
*
IMPLICIT NONE
*
* VARIABLES DEFINED IN UEL, PASSED BACK TO ABAQUS
*
REAL(8) :: RHS,AMATRX,SVARS,ENERGY
*
* VARIABLES PASSED INTO UEL
*
REAL(8) :: PROPS,coords,Uall,DUall,Vel,Accn,TIME,
1 DTIME,PARAMS,ADLMAG,PREDEF,DDLMAG,PNEWDT,PERIOD
INTEGER :: NDOFEL,NRHS,NSVARS,NPROPS,MCRD,NNODE,JTYPE,KSTEP,KINC,
1 JELEM,NDLOAD,JDLTYP,NPREDF,LFLAGS,MLVARX,MDLOAD,JPROPS,NJPROP
*
DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),PROPS(*),
1 SVARS(*),ENERGY(8),coords(MCRD,NNODE),Uall(NDOFEL),
2 DUall(MLVARX,*),Vel(NDOFEL),Accn(NDOFEL),TIME(2),PARAMS(*),
3 JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),
4 PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*)
if(jtype.eq.2) then
!
! This is a 2-D plane strain analysis, two nodes E-B beam element
!
nDim = 2
call UPE2(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,
+ PROPS,NPROPS,coords,MCRD,NNODE,Uall,DUall,Vel,Accn,JTYPE,
+ TIME,DTIME,KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,
+ PREDEF,NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,
+ NJPROP,PERIOD,nDim)
endif
return
end subroutine uel
*********
subroutine UPE2(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS,
+ PROPS,NPROPS,coords,MCRD,NNODE,Uall,DUall,Vel,Accn,JTYPE,
+ TIME,DTIME,KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,
+ PREDEF,NPREDF,LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,
+ NJPROP,PERIOD,nDim)
*
IMPLICIT NONE
*
* VARIABLES DEFINED IN UEL, PASSED BACK TO ABAQUS
*
REAL(8) :: RHS,AMATRX,SVARS,ENERGY
*
* VARIABLES PASSED INTO UEL
*
REAL(8) :: PROPS,coords,Uall,DUall,Vel,Accn,TIME,
1 DTIME,PARAMS,ADLMAG,PREDEF,DDLMAG,PNEWDT,PERIOD
INTEGER :: NDOFEL,NRHS,NSVARS,NPROPS,MCRD,NNODE,JTYPE,KSTEP,KINC,
1 JELEM,NDLOAD,JDLTYP,NPREDF,LFLAGS,MLVARX,MDLOAD,JPROPS,NJPROP,
1 ndim
*
DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),PROPS(*),
1 SVARS(*),ENERGY(8),coords(MCRD,NNODE),Uall(NDOFEL),
2 DUall(MLVARX,*),Vel(NDOFEL),Accn(NDOFEL),TIME(2),PARAMS(*),
3 JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*),
4 PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*)
real*8 vmass,wmass,susk,susc,rhob,areab,Eib,raydamp_a,raydamp_b,
+ x0,velo,vlen,Eab,GJb,eln,x1x,x2x,wgh,
+ xn,xc,xmax,xmin,stf(ndofel,ndofel),xksi,shph(ndofel,1),
+ phuu,phwu,mslp(ndofel,ndofel),a0,a1,a2,a3,a4,a5,a6,a7,
+ alpha,beta,gama,rdamp(ndofel,ndofel),ydx(ndofel),
+ gcrt(ndofel),glst(ndofel),fxt(ndofel),cntf,wdisp,wvel,wacc,
+ vinc,vfc,velt,qucrt,qwcrt,qulst,qwlst,qcrt(ndofel),
+ qlst(ndofel)
real*8 shcn(1,ndofel), shcn1d(1,ndofel),
* shcn2d(1,ndofel), bps, bps1d, bps2d
integer*4 nc,i,iflag,j,iter
real*8 zero,one,two,half,Pi,three,third,gacc,gauss_1d
real*8 chi(NNODE), aold
parameter(zero=0.d0,one=1.d0,two=2.d0,half=0.5d0,Pi=3.141592653d0,
+ three=3.d0,third=1.d0/3.d0,gacc=-9.81d0,gauss_1d=1.d0/sqrt(3.d0))
! Obtain material properties
rhob=props(1)
areab=props(2)
Eib=props(3)
Eab=props(4)
GJb=props(5)
vmass=props(6)
wmass=props(7)
susk=props(8)
susc=props(9)
x0=props(10)
velo=props(11)
vlen=props(12)
raydamp_a=props(13)
raydamp_b=props(14)
nc=jprops(1)
! Obtain time integration parameters, used in condensing the vehicle
! dofs into those of the wheel basing on the Newmark scheme
alpha=params(1)
beta=params(2)
gama=params(3)
a0=one/beta/dtime/dtime
a1=one/beta/dtime
a2=one/two/beta-one
a3=dtime*(one-gama)
a4=gama*dtime
a5=a4*a0
a6=a4*a1-one
a7=a4*a2-a3
! Suppose the vehicle is moving along x axle
! suppose spacing between neighboring vehicles are larger than
! the length of beam element, so that a given beam element will
! be loaded by ONLY one vehicle.
x1x=dabs(coords(1,1))
x2x=dabs(coords(1,2))
xmax=max(x1x,x2x)
xmin=min(x1x,x2x)
eln=xmax-xmin
if (eln .le. 0.0d0) then
write(80,*) 'Zero element length, pls check ...'
call xit
endif
iflag=0
iter=0
do i=1,nc,1
xn=x0-(i-1)*vlen
xc=xn+velo*time(2)
if (abs(xc).le.xmax .and. abs(xc).ge.xmin) then
! This conventional beam element will be transformed
! into VBI element
iter=i
iflag=1
xksi=(abs(xc)-xmin)/(xmax-xmin)
!write(80,*) xc,xmin,xksi
exit
endif
enddo
! Form stiffness matrix, upper half
stf=zero
stf(1,1)=Eab/eln
stf(1,4)=-Eab/eln
stf(2,2)=12.0d0*Eib/(eln**3)
stf(2,3)=6.0d0*Eib/(eln**2)
stf(2,5)=-12.0d0*Eib/(eln**3)
stf(2,6)=6.0d0*Eib/(eln**2)
stf(3,3)=4.0d0*Eib/eln
stf(3,5)=-6.0d0*Eib/(eln**2)
stf(3,6)=2.0d0*Eib/eln
stf(4,4)=Eab/eln
stf(5,5)=12.0d0*Eib/(eln**3)
stf(5,6)=-6.0d0*Eib/(eln**2)
stf(6,6)=4.0d0*Eib/eln
! Form mass matrix, lumped mass is used, upper half
!mslp=zero
!mslp(1,1)=rhob*areab*eln/two
!mslp(2,2)=rhob*areab*eln/two
!mslp(4,4)=rhob*areab*eln/two
!mslp(5,5)=rhob*areab*eln/two
! Form mass matrix, consist mass is used, upper half
wgh=rhob*areab*eln
mslp=zero
mslp(1,1)=third*wgh
mslp(1,4)=third*half*wgh
mslp(2,2)=156.0d0*wgh/420.0d0
mslp(2,3)=-22.0d0*eln*wgh/420.0d0
mslp(2,5)=54.0d0*wgh/420.0d0
mslp(2,6)=13.0d0*eln*wgh/420.0d0
mslp(3,3)=4.0d0*(eln**2)*wgh/420.0d0
mslp(3,5)=-13.0d0*eln*wgh/420.0d0
mslp(3,6)=-three*(eln**2)*wgh/420.0d0
mslp(4,4)=third*wgh
mslp(5,5)=156.0d0*wgh/420.0d0
mslp(5,6)=22.0d0*eln*wgh/420.0d0
mslp(6,6)=4.0d0*(eln**2)*wgh/420.0d0
! Form stiffness and mass matrix, lower half, use symmetry
do i=1,6,1
do j=i+1,6,1
stf(j,i)=stf(i,j)
mslp(j,i)=mslp(i,j)
enddo
enddo
! Form damp matrix, Rayleigh damping is assumed.
rdamp=raydamp_a*mslp+raydamp_b*stf
if (iflag .eq. 1) then
! Hermite interpolation for determing contribution
! of moving vehicle on beam stiffness matrix
! Hermite interpolation func at current local coord xksi
shph=zero
shph(2,1)=one-three*xksi**2+two*xksi**3
shph(5,1)=three*xksi**2-two*xksi**3
! Consistent interpolation
shcn=zero
shcn(1,2)=shph(2,1)
shcn(1,5)=shph(5,1)
shcn(1,3)=xksi*(one-two*xksi+xksi**2)*eln
shcn(1,6)=xksi*(xksi**2-xksi)*eln
! First order derivative of shcn with respect to xksi
shcn1d=zero
shcn1d(1,2)=two*three*xksi*(xksi-one)
shcn1d(1,5)=two*three*xksi*(one-xksi)
shcn1d(1,3)=(one-two*two*xksi+three*xksi**2)*eln
shcn1d(1,6)=(three*xksi**2-two*xksi)*eln
! Second order derivative of shcn with respect to xksi
shcn2d=zero
shcn2d(1,2)=two*three*(two*xksi-one)
shcn2d(1,5)=two*three*(one-two*xksi)
shcn2d(1,3)=(three*xksi-two)*two*eln
shcn2d(1,6)=(three*xksi-one)*two*eln
phuu=a0*vmass+a5*susc+susk
phwu=-a5*susc-susk
! Assembling stiffness for VBI element
! phwu, phuu should technically be matrices
!stf=stf+(one-phwu/phuu)*susk*matmul(shph,shcn)
! What are those?
!stf=stf+(wmass*(velo/eln)**2)*matmul(shph,shcn2d)
!stf=stf+(one+phwu/phuu)*susc*(velo/eln)*matmul(shph,shcn1d)
!rdamp=rdamp+(one-phwu/phuu)*susc*matmul(shph,shcn)
!rdamp=rdamp+(two*wmass*velo/eln)*matmul(shph,shcn1d)
rdamp = 0
! Update mass matrix for VBI element
! Just first term survives -> (mww): no triangular terms
!mslp=mslp+wmass*matmul(shph,shcn)
endif
if ((LFLAGS(1) .eq. 1) .or. (LFLAGS(1) .eq. 2)) then
if (LFLAGS(3) .eq. 1) then
AMATRX = stf
endif
endif
if (LFLAGS(1) .eq. 41 .or. LFLAGS(1) .eq. 99) then
if (LFLAGS(3) .eq. 2) then
AMATRX = stf
endif
if (LFLAGS(3) .eq. 3) then
AMATRX = rdamp
endif
if (LFLAGS(3) .eq. 4) then
AMATRX = mslp
endif
endif
! Time stepping scheme
if ((LFLAGS(1) .eq. 11) .or. (LFLAGS(1) .eq. 12)) then
AMATRX=a0*mslp+(one+alpha)*a5*rdamp+(one+alpha)*stf
endif
! SVARS(NSVARS) contains contribution of beam disp, beam velocity
! external force on beam, and contribution of moving vehicle only if
! the current beam element is a VBI element, at time t+deltat and t
! for internal variables at time t
! i.e. SVARS(1:6)=YDX|t=-C*Vel-K*Uall,svars(7:12)=FXT|t;svars(13:18)=q|t;
! svars(19:21)={deltau,vel_u,acc_u}|t; svars(22)=flag|t;
! for internal variables at time t-deltat
! i.e. SVARS(22+1:22+6)=YDX|t-deltat,svars(22+7:22+12)=FXT|t-deltat; &
! svars(22+13:22+18)=q|t-deltat,svars(22+19:22+21)={du,duvel,duacc}|t-deltat; &
! svars(22+22)=flag|t-deltat
! i.e. SVARS(NSVARS=44)
! G|t+deltat=-C*Vel|t+deltat-K*Uall|t+deltta+Fxt|t+deltat-svars(13:18)
! G|t=svars(1:6)+svars(7:12)-svars(35:40)
! RHS|t+deltat=-M*Acc|t+deltat+(1+alpha)*G|t+deltat-alpha*G|t
gcrt=zero
glst=zero
fxt=zero
ydx=zero
ydx=-matmul(rdamp,Vel)-matmul(stf,Uall)
glst=svars(1:6)
gcrt=ydx
svars(22)=real(iflag)
! Update just if at current iteration EL is VBI, otherwise {qc}=0
if (iflag .eq. 1) then
!xc - coordinate of the wheel at global coord sys
! bps: amplitude of road roughness at xc
! bps1d: the first derivative of road roughness at xc
! bps2d: the second derivative of road roughness at xc
xc=x0-(iter-1)*vlen+velo*time(2)
call roadrough (xc,bps,bps1d,bps2d)
do i=1,ndofel,1
fxt(i)=shph(i,1)*(vmass+wmass)*gacc
fxt(i)=fxt(i)-shph(i,1)*(one+phwu/phuu)*susk*bps
fxt(i)=fxt(i)-shph(i,1)*(one+phwu/phuu)*susc*bps1d*velo
fxt(i)=fxt(i)-shph(i,1)*wmass*bps2d*velo*velo
enddo
gcrt=gcrt+fxt
if (int(svars(22)) .eq. 1) then
gcrt=gcrt-svars(13:18)
glst=glst+svars(7:12)
endif
! if it also was at the previous increment, then q|t-deltat
! otherwise, q|t-deltat=0
if (int(svars(44)) .eq. 1) then
glst=glst-svars(35:40)
endif
! gcrt, glst have been calculated, so i can update svars
svars(29:43)=svars(7:21)
svars(7:12)=fxt
! Obtain disp, vel and acc at contact point for t+deltat
wdisp=zero
wvel=zero
wacc=zero
do i=1,ndofel,1
wdisp=wdisp+shcn(1,i)*Uall(i)
wvel=wvel+shcn(1,i)*Vel(i)
wacc=wacc+shcn(1,i)*Accn(i)
enddo
! Include contribution of road roughness into wheel response
!wdisp=wdisp+bps
!wvel=wvel+velo*bps1d
!wacc=wacc+velo*velo*bps2d
cntf=-susc*wvel-susk*wdisp
! forces on vehicle
vfc=vmass*(a1*svars(20)+a2*svars(21))+
+ susc*(a6*svars(20)+a7*svars(21))-susk*svars(19)
! eq.11
vinc=(vfc-cntf)/phuu
svars(19)=vinc
velt=svars(20)
aold=svars(21)
svars(21)=a0*vinc-a1*velt-a2*svars(21)
svars(20)=velt+a4*svars(21)+a3*aold
qucrt=vmass*(a1*svars(20)+a2*svars(21))+
+ susc*(a6*svars(20)+a7*svars(21))-susk*svars(19)
qwcrt=-susc*(a6*svars(20)+a7*svars(21))+susk*svars(19)
do i=1,ndofel,1
qcrt(i)=shph(i,1)*(phwu/phuu*qucrt-qwcrt)
enddo
svars(13:18)=qcrt
endif
!! Modify calculation in case of half step
if ((lflags(1) .eq. 11) .or. (lflags(1) .eq. 12)) then
rhs(1:ndofel,1)=-matmul(mslp,Accn)+(one+alpha)*gcrt-alpha*glst
endif
! Update internal variable at the end of current increment
svars(44)=svars(22)
svars(23:28)=svars(1:6)
svars(1:6)=ydx
! Update energy, energy(1): kinematic energy,
! energy(2): elastic strain energy
energy(1)=zero
do i=1,2,1
energy(1)=energy(1)+half*Vel(i)*(half*rhob*areab*eln)*Vel(i)
enddo
do i=4,5,1
energy(1)=energy(1)+half*Vel(i)*(half*rhob*areab*eln)*Vel(i)
enddo
energy(2)=zero
energy(2)=energy(2)+half*Eab*(Uall(1)-Uall(4))**2/eln
energy(2)=energy(2)+half*Eib*12.0d0*(Uall(2)-Uall(5))**2/eln**3
!Compute BM through Hermite interpolation
chi(1) = compute_curvature(-gauss_1d, eln, Uall)
chi(2) = compute_curvature(gauss_1d, eln, Uall)
svars(45)=-EIb*chi(1)
svars(46)=-EIb*chi(2)
return
contains
function compute_curvature(xi_gauss, L, d) result(curv)
implicit none
real(8), intent(in) :: xi_gauss ! Normalized position along the element (-1 <= xi <= 1)
real(8), intent(in) :: L ! Element length
real(8), intent(in) :: d(6) ! Nodal DOFs: [u1, w1, θ1, u2, w2, θ2]
real(8) :: curv ! Output: curvature w''(x)
real(8) :: xi ! Normalized position (0 <= xi <=1)
real(8) :: N1dd, N2dd, N3dd, N4dd
! get equivalent point in (0,1)
xi = (xi_gauss+1)/2
! Second derivatives of Hermite shape functions
N1dd = -6.0d0 + 12.0d0 * xi
N2dd = L * (-4.0d0 + 6.0d0 * xi)
N3dd = 6.0d0 - 12.0d0 * xi
N4dd = L * (-2.0d0 + 6.0d0 * xi)
! Compute curvature
curv = (1.0d0 / (L * L)) * (N1dd * d(2) + N2dd * d(3) +
+ N3dd * d(5) + N4dd * d(6))
end function compute_curvature
end subroutine UPE2
