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

Go to the source code of this file.

Functions/Subroutines

subroutine springstiff_n2f (xl, elas, konl, 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, ielorien, orab, norien, nelem)
 

Function/Subroutine Documentation

◆ springstiff_n2f()

subroutine springstiff_n2f ( real*8, dimension(3,10), intent(in)  xl,
real*8, dimension(21), intent(inout)  elas,
integer, dimension(20), intent(in)  konl,
real*8, dimension(0:mi(2),10), intent(in)  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(inout)  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, dimension(mi(3),*), intent(in)  ielorien,
real*8, dimension(7,*), intent(in)  orab,
integer, intent(in)  norien,
integer, intent(in)  nelem 
)
24 !
25 ! calculates the stiffness of a spring (node-to-face penalty)
26 !
27  implicit none
28 !
29  character*8 lakonl
30 !
31  integer konl(20),i,j,imat,ncmat_,ntmat_,k,l,nope,nterms,iflag,
32  & i1,kode,niso,id,nplicon(0:ntmat_,*),npmat_,nelcon(2,*),
33  & iperturb(*),nmethod,mi(*),ne0,nstate_,nasym,ielorien(mi(3),*),
34  & norien,nelem,idof,idof1,idof2,iorien
35 !
36  real*8 xl(3,10),elas(21),ratio(9),pproj(3),val,shp2(7,9),
37  & al(3),s(60,60),voldl(0:mi(2),10),pl(3,10),xn(3),dm,
38  & c1,c2,c3,c4,alpha,beta,elcon(0:ncmat_,ntmat_,*),xm(3),
39  & xmu(3,3,10),dxmu(3,10),dval(3,10),fpu(3,3,10),xi,et,
40  & xs2(3,7),t1l,elconloc(21),plconloc(802),xk,fk,
41  & xiso(200),yiso(200),dd0,plicon(0:2*npmat_,ntmat_,*),
42  & a11,a12,a22,b1(3,10),b2(3,10),dal(3,3,10),qxxy(3),fnl(3),
43  & qxyy(3),dxi(3,10),det(3,10),determinant,c11,c12,c22,
44  & qxyx(3),qyxy(3),springarea(2),dd,dist,t(3),tu(3,3,10),
45  & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),
46  & um,eps,pi,dftdt(3,3),tp(3),te(3),ftrial(3),clear,
47  & dftrial,dfnl,dfshear,dg,dte,alnew(3),dfn(3,10),reltime,
48  & overlap,pres,dpresdoverlap,overclosure,orab(7,*),
49  & xn1(3),xn2(3),a(3,3)
50 !
51  intent(in) xl,konl,voldl,imat,elcon,nelcon,ielorien,orab,
52  & ncmat_,ntmat_,nope,lakonl,t1l,kode,plicon,norien,nelem,
53  & nplicon,npmat_,iperturb,nmethod,mi,ne0,
54  & nstate_,xstateini,reltime,nasym
55 !
56  intent(inout) s,xstate,elas,springarea,elconloc
57 !
58  iflag=4
59 !
60 ! actual positions of the nodes belonging to the contact spring
61 ! (otherwise no contact force)
62 !
63  do i=1,nope
64  do j=1,3
65  pl(j,i)=xl(j,i)+voldl(j,i)
66  enddo
67  enddo
68 !
69  if(lakonl(7:7).eq.'A') then
70  if(lakonl(4:6).eq.'RNG') then
71 !
72 ! SPRINGA-element
73 !
74  dd0=dsqrt((xl(1,2)-xl(1,1))**2
75  & +(xl(2,2)-xl(2,1))**2
76  & +(xl(3,2)-xl(3,1))**2)
77  dd=dsqrt((pl(1,2)-pl(1,1))**2
78  & +(pl(2,2)-pl(2,1))**2
79  & +(pl(3,2)-pl(3,1))**2)
80  if(dd.le.0.d0) then
81  write(*,*) '*ERROR in springstiff_n2f: spring connects'
82  write(*,*) ' two nodes at the same location:'
83  write(*,*) ' spring length is zero'
84  call exit(201)
85  endif
86  do i=1,3
87  xn(i)=(pl(i,2)-pl(i,1))/dd
88  enddo
89  val=dd-dd0
90 !
91 ! interpolating the material data
92 !
93  call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l,
94  & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_)
95 !
96 ! calculating the spring force and the spring constant
97 !
98  if(kode.eq.2)then
99  xk=elconloc(1)
100  fk=xk*val
101  else
102  niso=int(plconloc(801))
103  do i=1,niso
104  xiso(i)=plconloc(2*i-1)
105  yiso(i)=plconloc(2*i)
106  enddo
107  call ident(xiso,val,niso,id)
108  if(id.eq.0) then
109  xk=0.d0
110  fk=yiso(1)
111  elseif(id.eq.niso) then
112  xk=0.d0
113  fk=yiso(niso)
114  else
115  xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
116  fk=yiso(id)+xk*(val-xiso(id))
117  endif
118  endif
119 !
120  c1=fk/dd
121  c2=xk-c1
122  do i=1,3
123  do j=1,3
124  s(i,j)=c2*xn(i)*xn(j)
125  enddo
126  s(i,i)=s(i,i)+c1
127  enddo
128  else
129 !
130 ! GAP-element
131 ! interpolating the material data
132 !
133  call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l,
134  & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_)
135 !
136  dd=elconloc(1)
137  xn(1)=elconloc(2)
138  xn(2)=elconloc(3)
139  xn(3)=elconloc(4)
140  xk=elconloc(5)
141  pi=4.d0*datan(1.d0)
142  eps=pi*elconloc(6)/xk
143  overclosure=-dd-xn(1)*(voldl(1,2)-voldl(1,1))
144  & -xn(2)*(voldl(2,2)-voldl(2,1))
145  & -xn(3)*(voldl(3,2)-voldl(3,1))
146  fk=-xk*overclosure*(0.5d0+datan(overclosure/eps)/pi)
147  c2=xk*((0.5d0+datan(overclosure/eps)/pi)
148  & +overclosure/(pi*eps*(1.d0+(overclosure/eps)**2)))
149  do i=1,3
150  do j=1,3
151  s(i,j)=c2*xn(i)*xn(j)
152  enddo
153  enddo
154  endif
155 !
156  do i=1,3
157  do j=1,3
158  s(i+3,j)=-s(i,j)
159  s(i,j+3)=-s(i,j)
160  s(i+3,j+3)=s(i,j)
161  enddo
162  enddo
163  return
164  elseif(lakonl(7:7).eq.'1') then
165 !
166 ! spring1-element
167 ! determine the direction of action of the spring
168 !
169  idof=nint(elcon(3,1,imat))
170 !
171  if(norien.gt.0) then
172  iorien=ielorien(1,nelem)
173  else
174  iorien=0
175  endif
176 !
177  if(iorien.eq.0) then
178  do i=1,3
179  xn(i)=0.d0
180  enddo
181  xn(idof)=1.d0
182  else
183  call transformatrix(orab(1,iorien),xl(1,1),a)
184  do i=1,3
185  xn(i)=a(i,idof)
186  enddo
187  endif
188 !
189 ! interpolating the material data
190 !
191  call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l,
192  & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_)
193 !
194 ! calculating the spring constant
195 !
196  if(kode.eq.2)then
197  xk=elconloc(1)
198  else
199  niso=int(plconloc(801))
200  do i=1,niso
201  xiso(i)=plconloc(2*i-1)
202  yiso(i)=plconloc(2*i)
203  enddo
204  call ident(xiso,val,niso,id)
205  if(id.eq.0) then
206  xk=0.d0
207  elseif(id.eq.niso) then
208  xk=0.d0
209  else
210  xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
211  endif
212  endif
213 !
214 ! assembling the stiffness matrix
215 !
216  do i=1,3
217  do j=1,3
218  s(i,j)=xk*xn(i)*xn(j)
219  enddo
220  enddo
221  return
222  elseif(lakonl(7:7).eq.'2') then
223 !
224 ! spring2-element
225 ! determine the direction of action of the spring
226 !
227  idof1=nint(elcon(3,1,imat))
228  idof2=nint(elcon(4,1,imat))
229 !
230  if(norien.gt.0) then
231  iorien=ielorien(1,nelem)
232  else
233  iorien=0
234  endif
235 !
236  if(iorien.eq.0) then
237  do i=1,3
238  xn1(i)=0.d0
239  xn2(i)=0.d0
240  enddo
241  xn1(idof1)=1.d0
242  xn2(idof2)=1.d0
243  else
244  call transformatrix(orab(1,iorien),xl(1,1),a)
245  do i=1,3
246  xn1(i)=a(i,idof1)
247  enddo
248  call transformatrix(orab(1,iorien),xl(1,2),a)
249  do i=1,3
250  xn2(i)=a(i,idof2)
251  enddo
252  endif
253 !
254 ! interpolating the material data
255 !
256  call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l,
257  & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_)
258 !
259 ! calculating the spring constant
260 !
261  if(kode.eq.2)then
262  xk=elconloc(1)
263  else
264  niso=int(plconloc(801))
265  do i=1,niso
266  xiso(i)=plconloc(2*i-1)
267  yiso(i)=plconloc(2*i)
268  enddo
269  call ident(xiso,val,niso,id)
270  if(id.eq.0) then
271  xk=0.d0
272  elseif(id.eq.niso) then
273  xk=0.d0
274  else
275  xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
276  endif
277  endif
278 !
279 ! assembling the stiffness matrix
280 !
281  do i=1,3
282  do j=1,3
283  s(i,j)=xk*xn1(i)*xn1(j)
284  s(i,j+3)=-xk*xn1(i)*xn2(j)
285  s(i+3,j)=-xk*xn2(i)*xn1(j)
286  s(i+3,j+3)=xk*xn2(i)*xn2(j)
287  enddo
288  enddo
289  return
290 !
291  endif
292 !
293 ! contact springs
294 !
295  nterms=nope-1
296 !
297 ! vector al connects the actual position of the slave node
298 ! with its projection on the master face
299 !
300  do i=1,3
301  pproj(i)=pl(i,nope)
302  enddo
303  call attach(pl,pproj,nterms,ratio,dist,xi,et)
304  do i=1,3
305  al(i)=pl(i,nope)-pproj(i)
306  enddo
307 !
308 ! determining the jacobian vector on the surface
309 !
310  if(nterms.eq.9) then
311  call shape9q(xi,et,pl,xm,xs2,shp2,iflag)
312  elseif(nterms.eq.8) then
313  call shape8q(xi,et,pl,xm,xs2,shp2,iflag)
314  elseif(nterms.eq.4) then
315  call shape4q(xi,et,pl,xm,xs2,shp2,iflag)
316  elseif(nterms.eq.6) then
317  call shape6tri(xi,et,pl,xm,xs2,shp2,iflag)
318  elseif(nterms.eq.7) then
319  call shape7tri(xi,et,pl,xm,xs2,shp2,iflag)
320  else
321  call shape3tri(xi,et,pl,xm,xs2,shp2,iflag)
322  endif
323 !
324 ! dxi(i,j) is the derivative of xi w.r.t. pl(i,j),
325 ! det(i,j) is the derivative of eta w.r.t. pl(i,j)
326 !
327 ! dxi and det are determined from the orthogonality
328 ! condition
329 !
330  a11=-(xs2(1,1)*xs2(1,1)+xs2(2,1)*xs2(2,1)+xs2(3,1)*xs2(3,1))
331  & +al(1)*xs2(1,5)+al(2)*xs2(2,5)+al(3)*xs2(3,5)
332  a12=-(xs2(1,1)*xs2(1,2)+xs2(2,1)*xs2(2,2)+xs2(3,1)*xs2(3,2))
333  & +al(1)*xs2(1,6)+al(2)*xs2(2,6)+al(3)*xs2(3,6)
334  a22=-(xs2(1,2)*xs2(1,2)+xs2(2,2)*xs2(2,2)+xs2(3,2)*xs2(3,2))
335  & +al(1)*xs2(1,7)+al(2)*xs2(2,7)+al(3)*xs2(3,7)
336 !
337  do i=1,3
338  do j=1,nterms
339  b1(i,j)=shp2(4,j)*xs2(i,1)-shp2(1,j)*al(i)
340  b2(i,j)=shp2(4,j)*xs2(i,2)-shp2(2,j)*al(i)
341  enddo
342  b1(i,nope)=-xs2(i,1)
343  b2(i,nope)=-xs2(i,2)
344  enddo
345 !
346  determinant=a11*a22-a12*a12
347  c11=a22/determinant
348  c12=-a12/determinant
349  c22=a11/determinant
350 !
351  do i=1,3
352  do j=1,nope
353  dxi(i,j)=c11*b1(i,j)+c12*b2(i,j)
354  det(i,j)=c12*b1(i,j)+c22*b2(i,j)
355  enddo
356  enddo
357 !
358 ! dal(i,j,k) is the derivative of al(i) w.r.t pl(j,k)
359 ! ( d al / d u_k)
360 !
361  do i=1,nope
362  do j=1,3
363  do k=1,3
364  dal(j,k,i)=-xs2(j,1)*dxi(k,i)-xs2(j,2)*det(k,i)
365  enddo
366  enddo
367  enddo
368  do i=1,nterms
369  do j=1,3
370  dal(j,j,i)=dal(j,j,i)-shp2(4,i)
371  enddo
372  enddo
373  do j=1,3
374  dal(j,j,nope)=dal(j,j,nope)+1.d0
375  enddo
376 !
377 ! d2q/dxx x dq/dy
378 !
379  qxxy(1)=xs2(2,5)*xs2(3,2)-xs2(3,5)*xs2(2,2)
380  qxxy(2)=xs2(3,5)*xs2(1,2)-xs2(1,5)*xs2(3,2)
381  qxxy(3)=xs2(1,5)*xs2(2,2)-xs2(2,5)*xs2(1,2)
382 !
383 ! dq/dx x d2q/dyy
384 !
385  qxyy(1)=xs2(2,1)*xs2(3,7)-xs2(3,1)*xs2(2,7)
386  qxyy(2)=xs2(3,1)*xs2(1,7)-xs2(1,1)*xs2(3,7)
387  qxyy(3)=xs2(1,1)*xs2(2,7)-xs2(2,1)*xs2(1,7)
388 !
389 ! Modified by Stefan Sicklinger
390 !
391 ! dq/dx x d2q/dxy
392 !
393  qxyx(1)=xs2(2,1)*xs2(3,6)-xs2(3,1)*xs2(2,6)
394  qxyx(2)=xs2(3,1)*xs2(1,6)-xs2(1,1)*xs2(3,6)
395  qxyx(3)=xs2(1,1)*xs2(2,6)-xs2(2,1)*xs2(1,6)
396 !
397 ! d2q/dxy x dq/dy
398 !
399  qyxy(1)=xs2(2,6)*xs2(3,2)-xs2(3,6)*xs2(2,2)
400  qyxy(2)=xs2(3,6)*xs2(1,2)-xs2(1,6)*xs2(3,2)
401  qyxy(3)=xs2(1,6)*xs2(2,2)-xs2(2,6)*xs2(1,2)
402 !
403 !
404 ! End modifications
405 !
406 ! normal on the surface
407 !
408  dm=dsqrt(xm(1)*xm(1)+xm(2)*xm(2)+xm(3)*xm(3))
409  do i=1,3
410  xn(i)=xm(i)/dm
411  enddo
412 !
413 ! distance from surface along normal (= clearance)
414 !
415  clear=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3)
416  if(nmethod.eq.1) then
417  clear=clear-springarea(2)*(1.d0-reltime)
418  endif
419 !
420 ! representative area: usually the slave surface stored in
421 ! springarea; however, if no area was assigned because the
422 ! node does not belong to any element, the master surface
423 ! is used
424 !
425  if(springarea(1).le.0.d0) then
426  if(nterms.eq.3) then
427  springarea(1)=dm/2.d0
428  else
429  springarea(1)=dm*4.d0
430  endif
431  endif
432 !
433 ! alpha and beta, taking the representative area into account
434 ! (conversion of pressure into force)
435 !
436  if(int(elcon(3,1,imat)).eq.1) then
437 !
438 ! exponential overclosure
439 !
440  if(dabs(elcon(2,1,imat)).lt.1.d-30) then
441  elas(1)=0.d0
442  elas(2)=0.d0
443  else
444  alpha=elcon(2,1,imat)*springarea(1)
445  beta=elcon(1,1,imat)
446  if(-beta*clear.gt.23.d0-dlog(alpha)) then
447  beta=(dlog(alpha)-23.d0)/clear
448  endif
449  elas(1)=dexp(-beta*clear+dlog(alpha))
450  elas(2)=-beta*elas(1)
451  endif
452  elseif(int(elcon(3,1,imat)).eq.2) then
453 !
454 ! linear overclosure
455 !
456  pi=4.d0*datan(1.d0)
457  eps=elcon(1,1,imat)*pi/elcon(2,1,imat)
458  elas(1)=(-springarea(1)*elcon(2,1,imat)*clear*
459  & (0.5d0+datan(-clear/eps)/pi))
460  elas(2)=-springarea(1)*elcon(2,1,imat)*
461  & ((0.5d0+datan(-clear/eps)/pi)-
462  & clear/(pi*eps*(1.d0+(clear/eps)**2)))
463  elseif(int(elcon(3,1,imat)).eq.3) then
464 !
465 ! tabular overclosure
466 !
467 ! interpolating the material data
468 !
469  call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l,
470  & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_)
471  overlap=-clear
472  niso=int(plconloc(801))
473  do i=1,niso
474  xiso(i)=plconloc(2*i-1)
475  yiso(i)=plconloc(2*i)
476  enddo
477  call ident(xiso,overlap,niso,id)
478  if(id.eq.0) then
479  dpresdoverlap=0.d0
480  pres=yiso(1)
481  elseif(id.eq.niso) then
482  dpresdoverlap=0.d0
483  pres=yiso(niso)
484  else
485  dpresdoverlap=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
486  pres=yiso(id)+dpresdoverlap*(overlap-xiso(id))
487  endif
488  elas(1)=springarea(1)*pres
489  elas(2)=-springarea(1)*dpresdoverlap
490  endif
491 !
492 ! contact force
493 !
494  do i=1,3
495  fnl(i)=-elas(1)*xn(i)
496  enddo
497 !
498 ! derivatives of the jacobian vector w.r.t. the displacement
499 ! vectors (d m / d u_k)
500 !
501  do k=1,nterms
502  xmu(1,1,k)=0.d0
503  xmu(2,2,k)=0.d0
504  xmu(3,3,k)=0.d0
505  xmu(1,2,k)=shp2(1,k)*xs2(3,2)-shp2(2,k)*xs2(3,1)
506  xmu(2,3,k)=shp2(1,k)*xs2(1,2)-shp2(2,k)*xs2(1,1)
507  xmu(3,1,k)=shp2(1,k)*xs2(2,2)-shp2(2,k)*xs2(2,1)
508  xmu(1,3,k)=-xmu(3,1,k)
509  xmu(2,1,k)=-xmu(1,2,k)
510  xmu(3,2,k)=-xmu(2,3,k)
511  enddo
512  do i=1,3
513  do j=1,3
514  xmu(i,j,nope)=0.d0
515  enddo
516  enddo
517 !
518 ! correction due to change of xi and eta
519 !
520  do k=1,nope
521  do j=1,3
522  do i=1,3
523 !
524 ! modified by Stefan Sicklinger
525 !
526  xmu(i,j,k)=xmu(i,j,k)+(qxxy(i)+qxyx(i))*dxi(j,k)
527  & +(qxyy(i)+qyxy(i))*det(j,k)
528  enddo
529  enddo
530  enddo
531 !
532 ! derivatives of the size of the jacobian vector w.r.t. the
533 ! displacement vectors (d ||m||/d u_k)
534 !
535  do k=1,nope
536  do i=1,3
537  dxmu(i,k)=xn(1)*xmu(1,i,k)+xn(2)*xmu(2,i,k)+
538  & xn(3)*xmu(3,i,k)
539  enddo
540 !
541 ! auxiliary variable: (d clear d u_k)*||m||
542 !
543  do i=1,3
544  dval(i,k)=al(1)*xmu(1,i,k)+al(2)*xmu(2,i,k)+
545  & al(3)*xmu(3,i,k)-clear*dxmu(i,k)+
546  & xm(1)*dal(1,i,k)+xm(2)*dal(2,i,k)+xm(3)*dal(3,i,k)
547  enddo
548 !
549  enddo
550 !
551  c1=1.d0/dm
552  c2=c1*c1
553  c3=elas(2)*c2
554  c4=elas(1)*c1
555 !
556 ! derivatives of the forces w.r.t. the displacement vectors
557 !
558  do k=1,nope
559  do j=1,3
560  do i=1,3
561  fpu(i,j,k)=-c3*xm(i)*dval(j,k)
562  & +c4*(xn(i)*dxmu(j,k)-xmu(i,j,k))
563  enddo
564  enddo
565  enddo
566 !
567 ! Coulomb friction for static calculations
568 !
569  if(ncmat_.ge.7) then
570  um=elcon(6,1,imat)
571  if(um.gt.0.d0) then
572 !
573 ! stiffness of shear stress versus slip curve
574 !
575  xk=elcon(7,1,imat)*springarea(1)
576 !
577 ! calculating the relative displacement between the slave node
578 ! and its projection on the master surface
579 !
580  do i=1,3
581  alnew(i)=voldl(i,nope)
582  do j=1,nterms
583  alnew(i)=alnew(i)-ratio(j)*voldl(i,j)
584  enddo
585  enddo
586 !
587 ! calculating the difference in relative displacement since
588 ! the start of the increment = lamda^*
589 !
590  do i=1,3
591  al(i)=alnew(i)-xstateini(3+i,1,ne0+konl(nope+1))
592  enddo
593 !
594 ! ||lambda^*||
595 !
596  val=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3)
597 !
598 ! update the relative tangential displacement
599 !
600  do i=1,3
601  t(i)=xstateini(6+i,1,ne0+konl(nope+1))+al(i)-val*xn(i)
602  enddo
603 !
604 ! store the actual relative displacement and
605 ! the actual relative tangential displacement
606 !
607  do i=1,3
608  xstate(3+i,1,ne0+konl(nope+1))=alnew(i)
609  xstate(6+i,1,ne0+konl(nope+1))=t(i)
610  enddo
611 !
612 ! d al/d u_k -> d al^*/d u_k
613 ! notice: xi & et are const.
614 !
615  do k=1,nope
616  do i=1,3
617  do j=1,3
618  dal(i,j,k)=0.d0
619  enddo
620  enddo
621  enddo
622 !
623  do i=1,nterms
624  do j=1,3
625  dal(j,j,i)=-shp2(4,i)
626  enddo
627  enddo
628 !
629  do j=1,3
630  dal(j,j,nope)=1.d0
631  enddo
632 !
633 ! (d al/d u_k).||m|| -> (d al^*/d u_k).||m||
634 !
635  do k=1,nope
636  do i=1,3
637  dval(i,k)=al(1)*xmu(1,i,k)+al(2)*xmu(2,i,k)
638  & +al(3)*xmu(3,i,k)-val*dxmu(i,k)
639  & +xm(1)*dal(1,i,k)+xm(2)*dal(2,i,k)
640  & +xm(3)*dal(3,i,k)
641  enddo
642  enddo
643 !
644 ! d t/d u_k
645 !
646  do k=1,nope
647  do j=1,3
648  do i=1,3
649  tu(i,j,k)=dal(i,j,k)
650  & -c1*(xn(i)*(dval(j,k)-val*dxmu(j,k))
651  & +val*xmu(i,j,k))
652  enddo
653  enddo
654  enddo
655 !
656 ! size of normal force
657 !
658  dfnl=dsqrt(fnl(1)**2+fnl(2)**2+fnl(3)**2)
659 !
660 ! maximum size of shear force
661 !
662  dfshear=um*dfnl
663 !
664 ! plastic and elastic slip
665 !
666  do i=1,3
667  tp(i)=xstateini(i,1,ne0+konl(nope+1))
668  te(i)=t(i)-tp(i)
669  enddo
670 !
671 ! the force due to normal contact is in -xn
672 ! direction (internal force) -> minus signs
673 !
674  do k=1,nope
675  do i=1,3
676  dfn(i,k)=-xn(1)*fpu(1,i,k)-xn(2)*fpu(2,i,k)-
677  & xn(3)*fpu(3,i,k)
678  enddo
679  enddo
680 !
681  dte=dsqrt(te(1)*te(1)+te(2)*te(2)+te(3)*te(3))
682 !
683 ! trial force
684 !
685  do i=1,3
686  ftrial(i)=xk*te(i)
687  enddo
688  dftrial=dsqrt(ftrial(1)**2+ftrial(2)**2+ftrial(3)**2)
689 !
690 ! check whether stick or slip
691 !
692  if((dftrial.lt.dfshear) .or. (dftrial.le.0.d0)) then
693 !
694 ! stick force
695 !
696  do i=1,3
697  fnl(i)=fnl(i)+ftrial(i)
698  xstate(i,1,ne0+konl(nope+1))=tp(i)
699  enddo
700 !
701 ! stick stiffness
702 !
703  do k=1,nope
704  do j=1,3
705  do i=1,3
706  fpu(i,j,k)=fpu(i,j,k)+xk*tu(i,j,k)
707  enddo
708  enddo
709  enddo
710  else
711 !
712 ! slip force
713 !
714  dg=(dftrial-dfshear)/xk
715  do i=1,3
716  ftrial(i)=te(i)/dte
717  fnl(i)=fnl(i)+dfshear*ftrial(i)
718  xstate(i,1,ne0+konl(nope+1))=tp(i)+dg*ftrial(i)
719  enddo
720 !
721 ! slip stiffness
722 !
723  c1=xk*dfshear/dftrial
724  do i=1,3
725  do j=1,3
726  dftdt(i,j)=-c1*ftrial(i)*ftrial(j)
727  enddo
728  dftdt(i,i)=dftdt(i,i)+c1
729  enddo
730 !
731  do k=1,nope
732  do j=1,3
733  do i=1,3
734  do l=1,3
735  fpu(i,j,k)=fpu(i,j,k)+dftdt(i,l)*tu(l,j,k)
736  enddo
737  if((nmethod.ne.4).or.(iperturb(1).gt.1)) then
738  fpu(i,j,k)=fpu(i,j,k)+um*ftrial(i)*dfn(j,k)
739  endif
740  enddo
741  enddo
742  enddo
743  endif
744  endif
745  endif
746 !
747 ! determining the stiffness matrix contributions
748 !
749 ! complete field shp2
750 !
751  shp2(1,nope)=0.d0
752  shp2(2,nope)=0.d0
753  shp2(4,nope)=-1.d0
754 !
755  do k=1,nope
756  do l=1,nope
757  do i=1,3
758  i1=i+(k-1)*3
759  do j=1,3
760  s(i1,j+(l-1)*3)=-shp2(4,k)*fpu(i,j,l)
761  & -shp2(1,k)*fnl(i)*dxi(j,l)
762  & -shp2(2,k)*fnl(i)*det(j,l)
763  enddo
764  enddo
765  enddo
766  enddo
767 !
768 ! symmetrizing the matrix
769 ! this is done in the absence of friction or for modal dynamic
770 ! calculations
771 !
772  if((nasym.eq.0).or.((nmethod.eq.4).and.(iperturb(1).le.1))) then
773  do j=1,3*nope
774  do i=1,j-1
775  s(i,j)=(s(i,j)+s(j,i))/2.d0
776  enddo
777  enddo
778  endif
779 !
780  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
dist
static double * dist
Definition: radflowload.c:42
transformatrix
subroutine transformatrix(xab, p, a)
Definition: transformatrix.f:20
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
a11
static double * a11
Definition: mafillkmain.c:30
b1
static double * b1
Definition: mafillkmain.c:30
shape6tri
subroutine shape6tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape6tri.f:20
attach
subroutine attach(pneigh, pnode, nterms, ratio, dist, xil, etl)
Definition: attach.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)