FEAP User Forum
FEAP => Programming => Topic started by: ad10 on July 09, 2021, 06:23:48 AM
-
Thanks for replying.
I am facing one issue here. I have formulated one 4-node simple linear elastic element. I am getting entirely different results for the following loop,
BATCh
DT,,1
LOOP,,10 ! Number of Time steps
TIME
!loop,,10 ! Number of NR iterations
Tang,,1
!NEXT
NEXT
ENd
If I don't use a time loop means simply Tang,,1, both results (FEAP and user element) are the same results. I can't be able to find the answer. Your suggestion is beneficial. Another question is related to the formulation of residual p,
I have seen two distinct cases for this,
1. by using stiffness
do i =1,8
p(i) = 0.d0
do j = 1,8
p(i) = p(i) - s(i,j)*ul(j)
end do
end do
For this I able to get results only for Tang,,1 Not for others (time loop)
2. by calculating, tr(B)*stress, but for this case, both look the same!
-
very hard to tell what the problem is.
All elements need to be formed as a residual in the p(*) array and a tangent in the s(*,*) array. A linear problem is just a special case.
If you include transients you have to use the ctan(3) values to pick up the time coefficient.
Post your element if you want a more specific answer.
-
Dear all,
I have formulated the S and p array in isw 3 as follows do i =1,ndim*nnode
do j = 1,ndim*nnode
s(i,j) = s(i,j) + kdd(i,j)*factor
end do
end do
do i =1,nst
do j = 1,nst
p(i) = p(i) - s(i,j)*utmp(j)
end do
end do
where, utmp are the current step displacements.
It works gives when only one NR iteration (simply Tang,,1) is used. But, it diverges for more NR loops used! What should I do so I will get a convergence?
-
Helpful commands:
SHOW,S -- tangent matrix in last element
SHOW,P -- RHS contribution from last element
SHOW,UL -- displacement vectors for last elementyou can use these command to check that you have the correct values in your arrays (check that they are consistent in matlab or python).
There is also the command
NTAN,ELEM,<element #>which shows you the difference between a numerical tangent and your element tangent (though I am not certain that this feature is working correctly in the current version of the code).
-
Post your elements we can look to see what you are doing that may cause the error
-
His element is attached to his last message: elmt03.f
-
Prof.,
I have checked for S and P.
In both cases, S is the same so I can say stiffness is calculations are right. But, p is different for both. In the case of displacement loading through edisp on edge, both results are the same! but results change for eforce/ force loading condition.
I have attached the S and P for both cases as follows.
-
Look at utmp values. They must be wrong?
-
Yeah, utmp values are wrong!
can you please tell me, what's I am missing here?
Thanks
-
How did you set utmp? Should come from ul array
-
Prof,
I have assigned the ul values to the utmp
! Initializing the displacement vector
utmp(1) = ul(1,1)
utmp(2) = ul(2,1)
! ---------------------------------
utmp(3) = ul(1,2)
utmp(4) = ul(2,2)
! ---------------------------------
utmp(5) = ul(1,3)
utmp(6) = ul(2,3)
! ---------------------------------
utmp(7) = ul(1,4)
utmp(8) = ul(2,4)
Like this and directly called these to formulate the load vector.
Is this the right way?
-
My isw 3 in short is this,
if(isw.eq.3 .or.isw.eq.6) then ! Compute residual/tangent
! Material Properties
emod = d(1)
enu = d(2)
thk = d(3)
nnode = nen
ndim = ndm
ngp = 2
! Initializing the displacement vector
utmp(1) = ul(1,1)
utmp(2) = ul(2,1)
! ---------------------------------
utmp(3) = ul(1,2)
utmp(4) = ul(2,2)
! ---------------------------------
utmp(5) = ul(1,3)
utmp(6) = ul(2,3)
! ---------------------------------
utmp(7) = ul(1,4)
utmp(8) = ul(2,4)
do i=1,nnode
do j=1,ndim
xref(j,i) = xl(j,i) !Reference coordinaates
end do
end do
! Initializing s(nst,nst) and p(nst)
do i =1,ndim*nnode
p(i) = zero
end do
do i =1,ndim*nnode
do j =1,ndim*nnode
s(i,j) = zero
end do
end do
! Setting the integration points
call int2d(ngp,lint,sg)
! Integration Loop
do l =1,lint
ss(1)=sg(1,l)
!write(*,*)'xi: ',ss(1)
ss(2)=sg(2,l)
!write(*,*)'eta: ',ss(2)
!write(*,*)'Step Integration: ',l
call kshape2d(shapef,dshape,ss(1),ss(2),ngp,nnode,ndim)
!write(*,*)'shapefunction: ',shapef
call kjac2d(xjac,xjacinv,detxjac,ndim,nnode,shapef
* ,dshape,xref)
call kstrain(stran,utmp,Bmat,ndim,nnode)
! Taking material model into consideration
call pzero(CmatVoigt2d,3*3)
CmatVoigt2D(1,1)=(emod/((one+enu)*(one-two*enu)))*
* (one-enu)
CmatVoigt2D(1,2)=(emod/((one+enu)*(one-two*enu)))*enu
CmatVoigt2D(2,1)=CmatVoigt2D(1,2)
CmatVoigt2D(2,2)=CmatVoigt2D(1,1)
CmatVoigt2D(3,3)=(emod/((one+enu)*(one-two*enu)))
* *((one-two*enu)/two)
call pzero(stress,3)
do ii = 1,3
do jj =1,3
stress(ii) = stress(ii) + CmatVoigt2d(ii,jj)*
* stran(jj)
end do
end do
!write(*,*)'CmatVoigt2D',CmatVoigt2D
factor = detxjac*sg(3,l)*thk
call kBmat2d(Bmat,ndim,nnode,ngp,shapef
* ,dshape,xjac,xjacinv)
!write(*,*)'Bmat: ',Bmat
call kCompTrans(Bmat,ndim*nnode,3,BmatT)
call kmatrixtripl(kdd,ndim*nnode,ndim*nnode,BmatT,
& nnode*ndim,3,CmatVoigt2d,3,3,Bmat,3,ndim*nnode)
!write(*,*)'kdd',kdd
!write(*,*)'kdd: ',kdd
do i =1,ndim*nnode
do j = 1,ndim*nnode
s(i,j) = s(i,j) + kdd(i,j)*factor
end do
end do
do i =1,nst
do j = 1,nst
p(i) = p(i) - s(i,j)*utmp(j)
end do
end do
end do
endif
I have nothing added to other ISWs. apart from 1.
-
Post the entire element for more help
-
I see it sorry
-
You have the computation of the p inside the quadrature loop, move it outside.
You do not need to initialize the p and s arrays, it is done before the call.
-
Thanks for pointing out, prof.
That's, why I am not getting the correct results. Thanks a lot! I will take note of it always.
-
You should study how feap stores elements in nodal blocks. Your element will not work if you have a problem where some element, say a shell, has 3 dofs/node.
Stiffness and residual are blocked in "ndf" sub-elements