Hi everybody,
i am trying to use feap 8.4 to implement a contact problem with friction. In this problem there are two workpieces, a upper workpiece and a lower workpiece as shown in initial.png in attachments. The lower workpiece is fixed. We move the upper workpiece first downwards to the lower workpiece and then towards the left.
If i use 3 nodes element, it works fine. You will find the inputfile Ineo3 and the plot neo3.png in attachments
But if i use 6 nodes element for upper workpiece and 3 nodes element for lower workpiece, you see the result is not good as shown in neo6.png. The inputfile is Ineo6.
As i know, feap does not support contact problem with higher order elements. So my idea for this problem with higher order element is that, first i implement Neo-Hooke material by myself as following, which works fine for no contact problem:
subroutine elmt37(d,ul,xl,ix,tl,s,r,ndf,ndm,nst,isw)
c
c Element with Neo-Hooke material
c
c 3-node Element / 1 Gauss-point
c property:
c
c lambda,mue for nonlinear material
c
c
implicit double precision (a-h,o-z)
c
include 'bdata.h'
include 'cdata.h' ! numnp,numel,nummat,nen,neq,ipr
include 'eldata.h'
include 'eltran.h'
include 'errchk.h'
include 'hdata.h'
include 'iofile.h'
include 'prstrs.h'
include 'comblk.h'
include 'tdata.h'
include 'strnum.h'
c
c
real*8 d(*),ul(ndf,nen),xl(ndm,*),tl(*),s(nst,*),r(*),
x shp(3,3),el(4,1),elg(3),wg,xsj,sig(3),eps(3)
integer ix(*)
real*8 lambda,mue,Bi(3,2),Bj(3,2),help(3,2),valp(13)
real*8 xd(2,3),kmat(2,2),kgeo(2,2),geo
real*8 gradu(2,2),c(3,3)
real*8 delta(2,2),sigt(2,2),xx,yy
real*8 defgrad(2,2),detf,detf2
real*8 lcgt(2,2)
integer i,j,k,i1,j1,k1,l,lint
c
c.... go to correct array processor
c 4 8
go to(1,2,3,3,2,3,2,3,2,2,2,2,2,3,3), isw
return
c.... input material properties
1 continue
c d(1) for Young's modulus
c d(2) for Poisson's ratio
call dinput(d(1),2)
lambda = d(1)*d(2)/((1.d0+d(2))*(1-2.d0*d(2)))
mue = d(1)/(2.d0*(1.d0+d(2)))
write(iow,2000) lambda,mue
istv = 13 ! Number of projection
return
2 continue
return
3 continue
c d(1) for Young's modulus
c d(2) for Poisson's ratio
lambda = d(1)*d(2)/((1.d0+d(2))*(1-2.d0*d(2)))
mue = d(1)/(2.d0*(1.d0+d(2)))
call tint2d(1,lint,el)
do 400 l = 1,lint
wg=el(4,l)
call trishp(el,xl,ndm,1,xsj,shp)
xsj=xsj*wg
xx=0.d0 !initial xx
yy=0.d0 !initial yy
do i=1,2
do j=1,2
gradu(i,j)=0.d0 !initial gradu
delta(i,j)=0.d0
if(i .eq. j) delta(i,j)=1.d0
enddo
enddo
do k=1,3
xx=xx+shp(3,k)*xl(1,k)
yy=yy+shp(3,k)*xl(2,k)
do i=1,2
xd(i,k)=xl(i,k)+ul(i,k) !xd for current position
do j=1,2
gradu(i,j) = gradu(i,j) + ul(i,k)*shp(j,k)
enddo
enddo
enddo
do i=1,2
do j=1,2
defgrad(i,j)=delta(i,j)+gradu(i,j) !deformation gradient
enddo
enddo
c determination of defgrad
detf=defgrad(1,1)*defgrad(2,2) - defgrad(1,2)*defgrad(2,1)
detf2=detf**2.d0
do i=1,2
do j=1,2
lcgt(i,j)=0.d0
do k=1,2
lcgt(i,j)=lcgt(i,j)+defgrad(i,k)*defgrad(j,k)
enddo
enddo
enddo
c left Cauchy-Green in eps
eps(1)=lcgt(1,1)
eps(2)=lcgt(2,2)
eps(3)=lcgt(1,2)
call hooke_act_37(eps,detf2,lambda,mue,sig,c)
sigt(1,1)=sig(1)
sigt(1,2)=sig(3)
sigt(2,1)=sig(3)
sigt(2,2)=sig(2)
if (mod(isw,3).eq.0) then
call trishp(el,xd,ndm,1,xsj,shp)
xsj=xsj*wg
do 300 i=1,3
Bi(1,1)=shp(1,i)
Bi(1,2)=0.d0
Bi(2,1)=0.d0
Bi(2,2)=shp(2,i)
Bi(3,1)=shp(2,i)
Bi(3,2)=shp(1,i)
if (isw.eq.3) then
do 310 j=1,3
Bj(1,1)=shp(1,j)
Bj(1,2)=0.d0
Bj(2,1)=0.d0
Bj(2,2)=shp(2,j)
Bj(3,1)=shp(2,j)
Bj(3,2)=shp(1,j)
geo=shp(1,i)*shp(1,j)*sig(1)
x +shp(2,i)*shp(2,j)*sig(2)
x +(shp(1,i)*shp(2,j)+shp(2,i)*shp(1,j))*sig(3)
do i1=1,3
do j1=1,2
help(i1,j1)=0.d0
do k1=1,3
help(i1,j1)=help(i1,j1)+c(i1,k1)*Bj(k1,j1)
enddo
enddo
enddo
do i1=1,2
do j1=1,2
kmat(i1,j1)=0.d0
kgeo(i1,j1)=delta(i1,j1)*geo*xsj
do k1=1,3
kmat(i1,j1)=kmat(i1,j1)+Bi(k1,i1)*help(k1,j1)*xsj
enddo
enddo
enddo
c compute K-term
do i1=1,2
do j1=1,2
s(2*(i-1)+i1,2*(j-1)+j1)=s(2*(i-1)+i1,2*(j-1)+j1)
x +kmat(i1,j1)+kgeo(i1,j1)
enddo
enddo
310 continue
endif
c compute residual -(Bt sigma)
do i1=1,2
do k1=1,3
r(2*(i-1)+i1)=r(2*(i-1)+i1)-Bi(k1,i1)*sig(k1)*xsj
enddo
enddo
300 continue
endif
if(isw.eq.4) then
write(iow,2500) l,xx,yy,defgrad,sig
write(*,2500) l,xx,yy,defgrad,sig
endif
if (isw.eq.8) then
valp(1)=defgrad(1,1)
valp(2)=defgrad(2,2)
valp(3)=defgrad(1,2)
valp(4)=defgrad(2,1)
valp(5)=sig(1)
valp(6)=sig(2)
valp(7)=sig(3)
valp(8)=0.d0
valp(9)=0.d0
valp(10)=0.d0
valp(11)=0.d0
valp(12)=0.d0
valp(13)=0.d0
call elmt37plot(ix,valp,shp,xsj,
x hr(nph),hr(nph+numnp),hr(ner),
x nel,numnp)
endif
400 continue
return
c Format
2000 format(' 3-Node-Element for Large Strain ',//,
x ' (Formulation in actual configuration) ',//,
x ' lambda = ',3e12.4,/,
x ' mue = ',3e12.4,/,
x ' ',3e12.4,/)
2500 format(' Element solution (GP-NR) : ',i4,/,
x ' x_1, x_2 (ref. conf.) : ',2e20.4,/,
x ' def1, def2, def3, def4 : ',4e20.4,/,
x ' sig_11, sig_22 , sig_12 : ',3e12.4,/,
x ' S_11, S_12, S_21, S_22 : ',4e12.4,/)
3000 format(' Warning: Element called with isw=',i4,' ! ',/,
x ' No action implemented ! ',/)
end
c-----7-----------------------------------------------------------------
subroutine elmt37plot(ix,value,shp,xsj,dt,st,ser,nel,numnp)
implicit none
integer nel,nen,numnp
integer i,l,ll
real*8 xg
integer ix(*)
real*8 dt(numnp),st(numnp,*),ser(*)
real*8 xsj,shp(3,3),value(13)
integer k
save
c Lumped and consistent projection routine
c Compute lumped projection and assemble stress integrals
do i = 1,nel
ll = ix(i)
if(ll.gt.0) then
xg = shp(3,i)*xsj
dt(ll) = dt(ll) + xg
c Projections of GP-Data
do k=1,13
st(ll,k) = st(ll,k) + value(k)*xg
enddo
endif
enddo
end
c-----7-----------------------------------------------------------------
subroutine hooke_act_37(eps,detf2,lambda,mue,sig,c)
c
c IN: eps - left Cauchy-Green-Tensor
c detf2 - (det F)^2
c lambda,mue - Materialconstanten
c OUT: sig - Cauchy-stress
c c - Materialtangent
c
real*8 eps(3),lambda,mue
real*8 sig(3),c(3,3),detf2
real*8 fakt1,fakt2,fakt3,detf
detf=dsqrt(detf2)
fakt1=0.5d0*lambda*(detf2-1.d0)/detf
fakt2=(lambda+2.d0*mue)/detf
fakt3=lambda*detf
sig(1)=mue/detf*(eps(1)-1.d0)+fakt1
sig(2)=mue/detf*(eps(2)-1.d0)+fakt1
sig(3)=mue/detf*eps(3)
c(1,1)=fakt2
c(1,2)=fakt3
c(1,3)=0.d0
c(2,1)=c(1,2)
c(2,2)=fakt2
c(2,3)=0.d0
c(3,1)=c(1,3)
c(3,2)=c(2,3)
c(3,3)=mue/detf-fakt1
end
c-----7-----------------------------------------------------------------
I choose node to segment strategy. Let the nodes of the upper workpiece be the slave nodes und the upper boundary of the lower workpiece be the master segment. So we can write the weak form for this problem and derive the tangent matrix and the residuum, correspondingly i will modify the element tangent matrix s and element residuum r in subroutine elmt37(d,ul,xl,ix,tl,s,r,ndf,ndm,nst,isw). Because the geometry of the master line is very simple and just a fixed straight line, so i think it is not a very hard task to implement. If i can finish this, we don't need to use the contact module of the feap, we solve this problem just like a normal problem without contact. How do you think my idea? Is NTOS a good strategy, or i should use segment to segment strategy for higher oder problem?
Best regards
