1!
2!     CalculiX - A 3-dimensional finite element program
3!              Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19      subroutine springstiff_f2f(xl,elas,voldl,s,imat,elcon,nelcon,
20     &  ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc,plicon,
21     &  nplicon,npmat_,iperturb,springarea,nmethod,mi,ne0,
22     &  nstate_,xstateini,xstate,reltime,nasym,
23     &  jfaces,igauss,pslavsurf,pmastsurf,clearini,kscale)
24!
25!     calculates the stiffness of a spring (face-to-face penalty)
26!
27      implicit none
28!
29      character*8 lakonl
30!
31      integer i,j,imat,ncmat_,ntmat_,k,l,nope,iflag,
32     &  kode,niso,id,nplicon(0:ntmat_,*),npmat_,nelcon(2,*),
33     &  iperturb(*),nmethod,mi(*),ne0,nstate_,nasym,
34     &  jfaces,igauss,nopem,nopes,nopep,kscale
35!
36      real*8 xl(3,19),elas(21),pproj(3),val,shp2m(7,9),
37     &  al(3),s(60,60),voldl(0:mi(2),19),pl(3,19),xn(3),
38     &  c1,c3,alpha,beta,elcon(0:ncmat_,ntmat_,*),xm(3),
39     &  fpu(3,3),xi,et,fnl(3),
40     &  xs2(3,7),t1l,elconloc(21),plconloc(82),xk,stickslope,
41     &  xiso(20),yiso(20),plicon(0:2*npmat_,ntmat_,*),
42     &  springarea(2),t(3),tu(3,3),overlap,pres,dpresdoverlap,
43     &  xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),
44     &  um,dftdt(3,3),tp(3),te(3),ftrial(3),clear,
45     &  dftrial,dfnl,dfshear,dg,dte,alnew(3),dfn(3),reltime,
46     &  xsj2s(3),xs2s(3,7),shp2s(7,9),weight,pslavsurf(3,*),
47     &  pmastsurf(6,*),clearini(3,9,*)
48!
49      iflag=1
50!
51!     # of master nodes
52!
53      nopem=ichar(lakonl(8:8))-48
54!
55!     # of slave nodes
56!
57      nopes=nope-nopem
58!
59!     actual positions of the nodes belonging to the contact spring
60!     (otherwise no contact force)
61!
62!     master nodes
63!
64      do i=1,nopem
65         do j=1,3
66            pl(j,i)=xl(j,i)+voldl(j,i)
67         enddo
68      enddo
69!
70!     slave nodes
71!
72      do i=nopem+1,nope
73         do j=1,3
74            pl(j,i)=xl(j,i)+voldl(j,i)+clearini(j,i-nopem,jfaces)
75     &             *reltime
76         enddo
77      enddo
78!
79!     contact springs
80!
81      xi=pslavsurf(1,igauss)
82      et=pslavsurf(2,igauss)
83      weight=pslavsurf(3,igauss)
84!
85      if(nopes.eq.8) then
86          call shape8q(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
87      elseif(nopes.eq.4) then
88          call shape4q(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
89      elseif(nopes.eq.6) then
90          call shape6tri(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
91      else
92          call shape3tri(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag)
93      endif
94!
95      nopep=nope+1
96!
97!     position and displacements of the integration point in the
98!     slave face
99!
100      do k=1,3
101          pl(k,nopep)=0.d0
102          voldl(k,nopep)=0.d0
103          do j=1,nopes
104              pl(k,nopep)=pl(k,nopep)+shp2s(4,j)*pl(k,nopem+j)
105              voldl(k,nopep)=voldl(k,nopep)+shp2s(4,j)*voldl(k,nopem+j)
106          enddo
107      enddo
108!
109      xi=pmastsurf(1,igauss)
110      et=pmastsurf(2,igauss)
111!
112!     determining the jacobian vector on the surface
113!
114      if(nopem.eq.8) then
115         call shape8q(xi,et,pl,xm,xs2,shp2m,iflag)
116      elseif(nopem.eq.4) then
117         call shape4q(xi,et,pl,xm,xs2,shp2m,iflag)
118      elseif(nopem.eq.6) then
119         call shape6tri(xi,et,pl,xm,xs2,shp2m,iflag)
120      else
121         call shape3tri(xi,et,pl,xm,xs2,shp2m,iflag)
122      endif
123!
124!     position of the projection of the slave integration point
125!     on the master faces (done at the start of the increment)
126!
127      do i=1,3
128         pproj(i)=0.d0
129         do j=1,nopem
130            pproj(i)=pproj(i)+shp2m(4,j)*pl(i,j)
131         enddo
132      enddo
133!
134!     vector connecting the integration point with its projection
135!     on the master face
136!
137      do i=1,3
138         al(i)=pl(i,nopep)-pproj(i)
139      enddo
140!
141!     normal vector on master face
142!
143      xn(1)=pmastsurf(4,igauss)
144      xn(2)=pmastsurf(5,igauss)
145      xn(3)=pmastsurf(6,igauss)
146!
147!     distance from surface along normal (= clearance)
148!
149      clear=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3)
150      if(nmethod.eq.1) then
151         clear=clear-springarea(2)*(1.d0-reltime)
152      endif
153!
154!     alpha and beta, taking the representative area into account
155!     (conversion of pressure into force)
156!
157      if(int(elcon(3,1,imat)).eq.1) then
158!
159!        exponential overclosure
160!
161         if(dabs(elcon(2,1,imat)).lt.1.d-30) then
162            elas(1)=0.d0
163            elas(2)=0.d0
164         else
165            alpha=elcon(2,1,imat)*springarea(1)
166            beta=elcon(1,1,imat)
167            if(-beta*clear.gt.23.d0-dlog(alpha)) then
168               beta=(dlog(alpha)-23.d0)/clear
169            endif
170            elas(1)=dexp(-beta*clear+dlog(alpha))
171            elas(2)=-beta*elas(1)
172         endif
173      elseif((int(elcon(3,1,imat)).eq.2).or.
174     &       (int(elcon(3,1,imat)).eq.4)) then
175!
176!        linear overclosure/tied overclosure
177!
178c        write(*,*) 'springstiff_f2f',springarea(1),elcon(2,1,imat),
179c     &       kscale
180         elas(2)=-springarea(1)*elcon(2,1,imat)/kscale
181         elas(1)=elas(2)*clear
182      elseif(int(elcon(3,1,imat)).eq.3) then
183!
184!        tabular overclosure
185!
186!        interpolating the material data
187!
188         call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l,
189     &     elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_)
190         overlap=-clear
191         niso=int(plconloc(81))
192         do i=1,niso
193            xiso(i)=plconloc(2*i-1)
194            yiso(i)=plconloc(2*i)
195         enddo
196         call ident(xiso,overlap,niso,id)
197         if(id.eq.0) then
198            dpresdoverlap=0.d0
199            pres=yiso(1)
200         elseif(id.eq.niso) then
201            dpresdoverlap=0.d0
202            pres=yiso(niso)
203         else
204            dpresdoverlap=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
205            pres=yiso(id)+dpresdoverlap*(overlap-xiso(id))
206         endif
207         elas(1)=springarea(1)*pres
208         elas(2)=-springarea(1)*dpresdoverlap
209      endif
210!
211!     contact force
212!
213      do i=1,3
214         fnl(i)=-elas(1)*xn(i)
215      enddo
216!
217      c3=elas(2)
218!
219!     derivatives of the forces w.r.t. the displacement vectors
220!
221      do j=1,3
222         do i=1,3
223            fpu(i,j)=-c3*xn(i)*xn(j)
224         enddo
225      enddo
226!
227!     Coulomb friction for static calculations
228!
229      if((ncmat_.ge.7).or.(int(elcon(3,1,imat)).eq.4)) then
230!
231!        tied contact
232!
233         if(int(elcon(3,1,imat)).eq.4) then
234            um=1.d30
235         else
236            um=elcon(6,1,imat)
237         endif
238         stickslope=elcon(7,1,imat)/kscale
239!
240         if(um.gt.0.d0) then
241!
242!     stiffness of shear stress versus slip curve
243!
244            xk=stickslope*springarea(1)
245!
246!     calculating the relative displacement between the slave node
247!     and its projection on the master surface
248!
249            do i=1,3
250               alnew(i)=voldl(i,nopep)
251               do j=1,nopem
252                  alnew(i)=alnew(i)-shp2m(4,j)*voldl(i,j)
253               enddo
254            enddo
255!
256!     calculating the difference in relative displacement since
257!     the start of the increment = lamda^*
258!
259            do i=1,3
260               al(i)=alnew(i)-xstateini(3+i,1,ne0+igauss)
261            enddo
262!
263!     ||lambda^*||
264!
265            val=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3)
266!
267!     update the relative tangential displacement
268!
269            do i=1,3
270               t(i)=xstateini(6+i,1,ne0+igauss)+al(i)-val*xn(i)
271            enddo
272!
273!     store the actual relative displacement and
274!     the actual relative tangential displacement
275!
276            do i=1,3
277               xstate(3+i,1,ne0+igauss)=alnew(i)
278               xstate(6+i,1,ne0+igauss)=t(i)
279            enddo
280!
281!     d t/d u_k
282!
283            do j=1,3
284               do i=1,3
285                  tu(i,j)=-xn(i)*xn(j)
286               enddo
287               tu(j,j)=tu(j,j)+1.d0
288            enddo
289!
290!     size of normal force
291!
292            dfnl=dsqrt(fnl(1)**2+fnl(2)**2+fnl(3)**2)
293!
294!     maximum size of shear force
295!
296            if(int(elcon(3,1,imat)).eq.4) then
297               dfshear=1.d30
298            else
299               dfshear=um*dfnl
300            endif
301c            dfshear=um*dfnl
302!
303!     plastic and elastic slip
304!
305            do i=1,3
306               tp(i)=xstateini(i,1,ne0+igauss)
307               te(i)=t(i)-tp(i)
308            enddo
309!
310            dte=dsqrt(te(1)*te(1)+te(2)*te(2)+te(3)*te(3))
311!
312!     trial force
313!
314            do i=1,3
315               ftrial(i)=xk*te(i)
316            enddo
317            dftrial=dsqrt(ftrial(1)**2+ftrial(2)**2+ftrial(3)**2)
318!
319!     check whether stick or slip
320!
321            if((dftrial.lt.dfshear) .or. (dftrial.le.0.d0)) then
322!
323!     stick force
324!
325               do i=1,3
326                  fnl(i)=fnl(i)+ftrial(i)
327                  xstate(i,1,ne0+igauss)=tp(i)
328               enddo
329!
330!     stick stiffness
331!
332               do j=1,3
333                  do i=1,3
334                     fpu(i,j)=fpu(i,j)+xk*tu(i,j)
335                  enddo
336               enddo
337            else
338!
339!     slip force
340!
341               dg=(dftrial-dfshear)/xk
342               do i=1,3
343                  ftrial(i)=te(i)/dte
344                  fnl(i)=fnl(i)+dfshear*ftrial(i)
345                  xstate(i,1,ne0+igauss)=tp(i)+dg*ftrial(i)
346               enddo
347!
348!     slip stiffness
349!
350               do i=1,3
351                  dfn(i)=-xn(1)*fpu(1,i)-xn(2)*fpu(2,i)-
352     &                 xn(3)*fpu(3,i)
353               enddo
354!
355               c1=xk*dfshear/dftrial
356               do i=1,3
357                  do j=1,3
358                     dftdt(i,j)=-c1*ftrial(i)*ftrial(j)
359                  enddo
360                  dftdt(i,i)=dftdt(i,i)+c1
361               enddo
362!
363               do j=1,3
364                  do i=1,3
365                     do l=1,3
366                        fpu(i,j)=fpu(i,j)+
367     &                                 dftdt(i,l)*tu(l,j)
368                     enddo
369                     if((nmethod.ne.4).or.(iperturb(1).gt.1)) then
370                        fpu(i,j)=fpu(i,j)+
371     &                                 um*ftrial(i)*dfn(j)
372                     endif
373                  enddo
374               enddo
375            endif
376         endif
377      endif
378!
379!     determining the stiffness matrix contributions
380!
381!     dFkm/dUlm
382!
383      do k=1,nopem
384         do i=1,3
385             do l=1,nopem
386               do j=1,3
387                  s(i+(k-1)*3,j+(l-1)*3)=
388     &                 shp2m(4,k)*shp2m(4,l)*fpu(i,j)
389               enddo
390            enddo
391         enddo
392      enddo
393!
394!     dFks/dUls
395!
396      do k=nopem+1,nopem+nopes
397         do i=1,3
398            do l=nopem+1,nopem+nopes
399               do j=1,3
400                  s(i+(k-1)*3,j+(l-1)*3)=shp2s(4,k-nopem)*
401     &                                shp2s(4,l-nopem)*fpu(i,j)
402               enddo
403            enddo
404         enddo
405      enddo
406!
407!     dFkm/dUls
408!
409      do k=1,nopem
410         do i=1,3
411            do l=nopem+1,nopem+nopes
412               do j=1,3
413                  s(i+(k-1)*3,j+(l-1)*3)=-shp2s(4,l-nopem)*
414     &                shp2m(4,k)*fpu(i,j)
415               enddo
416            enddo
417         enddo
418      enddo
419!
420!     dFks/dUlm
421!
422      do k=nopem+1,nopem+nopes
423         do i=1,3
424            do l=1,nopem
425               do j=1,3
426                  s(i+(k-1)*3,j+(l-1)*3)=
427     &                -shp2s(4,k-nopem)*shp2m(4,l)*fpu(i,j)
428               enddo
429            enddo
430         enddo
431      enddo
432!
433!     symmetrizing the matrix
434!     this is done in the absence of friction or for modal dynamic
435!     calculations
436!
437      if((nasym.eq.0).or.((nmethod.eq.4).and.(iperturb(1).le.1))) then
438         do j=1,3*nope
439            do i=1,j-1
440               s(i,j)=(s(i,j)+s(j,i))/2.d0
441            enddo
442         enddo
443      endif
444!
445      return
446      end
447
448