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