CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
springstiff_f2f.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine springstiff_f2f (xl, elas, voldl, s, imat, elcon, nelcon, ncmat_, ntmat_, nope, lakonl, t1l, kode, elconloc, plicon, nplicon, npmat_, iperturb, springarea, nmethod, mi, ne0, nstate_, xstateini, xstate, reltime, nasym, iloc, jfaces, igauss, pslavsurf, pmastsurf, clearini, kscale)
 

Function/Subroutine Documentation

◆ springstiff_f2f()

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