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 springforc_f2f(xl,vl,imat,elcon,nelcon, 20 & elas,fnl,ncmat_,ntmat_,nope,lakonl,t1l,kode,elconloc, 21 & plicon,nplicon,npmat_,senergy,nener,cstr,mi, 22 & springarea,nmethod,ne0,nstate_,xstateini, 23 & xstate,reltime,ielas,jfaces,igauss, 24 & pslavsurf,pmastsurf,clearini,venergy,kscale, 25 & konl,iout,nelem,smscale,mscalmethod) 26! 27! calculates the force of the spring (face-to-face penalty) 28! 29 implicit none 30! 31 character*8 lakonl 32! 33 integer i,j,k,imat,ncmat_,ntmat_,nope,iflag,mi(*), 34 & kode,niso,id,nplicon(0:ntmat_,*),npmat_,nelcon(2,*),nener, 35 & nmethod,ne0,nstate_,ielas,jfaces,kscale,konl(26), 36 & igauss,nopes,nopem,nopep,iout,nelem,mscalmethod 37! 38 real*8 xl(3,10),elas(21),t1l,al(3),vl(0:mi(2),19),stickslope, 39 & pl(3,19),xn(3),alpha,beta,fnl(3,19),tp(3),te(3),ftrial(3), 40 & t(3),dftrial,elcon(0:ncmat_,ntmat_,*),pproj(3),clear, 41 & xi,et,elconloc(*),plconloc(82),xk,val,xiso(20),yiso(20), 42 & plicon(0:2*npmat_,ntmat_,*),um,senergy,cstr(6),dg,venergy, 43 & dfshear,dfnl,springarea(2),overlap,pres,clearini(3,9,*), 44 & xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),t1(3),t2(3), 45 & dt1,dte,alnew(3),reltime,weight,xsj2m(3),xs2m(3,7),shp2m(7,9), 46 & xsj2s(3),xs2s(3,7),shp2s(7,9),pslavsurf(3,*),pmastsurf(6,*), 47 & smscale(*) 48! 49 include "gauss.f" 50! 51 iflag=2 52! 53 if(nener.eq.1) venergy=0.d0 54! 55! # of master nodes 56! 57 nopem=ichar(lakonl(8:8))-48 58! 59! # of slave nodes 60! 61 nopes=nope-nopem 62! 63! actual positions of the master nodes belonging to the contact spring 64! 65 do i=1,nopem 66 do j=1,3 67 pl(j,i)=xl(j,i)+vl(j,i) 68 enddo 69 enddo 70! 71! actual positions of the slave nodes belonging to the contact spring 72! 73 if(nmethod.ne.2) then 74 do i=nopem+1,nope 75 do j=1,3 76 pl(j,i)=xl(j,i)+clearini(j,i-nopem,jfaces)*reltime 77 & +vl(j,i) 78 enddo 79 enddo 80 else 81! 82! for frequency calculations the eigenmodes are freely 83! scalable, leading to problems with contact finding 84! 85 do i=nopem+1,nope 86 do j=1,3 87 pl(j,i)=xl(j,i)+clearini(j,i-nopem,jfaces)*reltime 88 enddo 89 enddo 90 endif 91! 92! location of integration point in slave face 93! 94 xi=pslavsurf(1,igauss) 95 et=pslavsurf(2,igauss) 96 weight=pslavsurf(3,igauss) 97! 98 iflag=1 99 if(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 else 106 call shape3tri(xi,et,pl(1,nopem+1),xsj2s,xs2s,shp2s,iflag) 107 endif 108! 109 nopep=nope+1 110! 111 do k=1,3 112 pl(k,nopep)=0.d0 113 vl(k,nopep)=0.d0 114 do j=1,nopes 115 pl(k,nopep)=pl(k,nopep)+shp2s(4,j)*pl(k,nopem+j) 116 vl(k,nopep)=vl(k,nopep)+shp2s(4,j)*vl(k,nopem+j) 117 enddo 118 enddo 119! 120! corresponding location in the master face 121! 122 xi=pmastsurf(1,igauss) 123 et=pmastsurf(2,igauss) 124! 125! determining the jacobian vector on the surface 126! 127 iflag=2 128 if(nopem.eq.8) then 129 call shape8q(xi,et,pl,xsj2m,xs2m,shp2m,iflag) 130 elseif(nopem.eq.4) then 131 call shape4q(xi,et,pl,xsj2m,xs2m,shp2m,iflag) 132 elseif(nopem.eq.6) then 133 call shape6tri(xi,et,pl,xsj2m,xs2m,shp2m,iflag) 134 else 135 call shape3tri(xi,et,pl,xsj2m,xs2m,shp2m,iflag) 136 endif 137! 138 do i=1,3 139 pproj(i)=0.d0 140 do j=1,nopem 141 pproj(i)=pproj(i)+shp2m(4,j)*pl(i,j) 142 enddo 143 enddo 144! 145! distance vector between both 146! 147 do i=1,3 148 al(i)=pl(i,nopep)-pproj(i) 149 enddo 150! 151! normal on the master face 152! 153 xn(1)=pmastsurf(4,igauss) 154 xn(2)=pmastsurf(5,igauss) 155 xn(3)=pmastsurf(6,igauss) 156! 157! distance from surface along normal (= clearance) 158! 159 clear=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3) 160! 161! check for a reduction of the initial penetration, if any 162! 163 if(nmethod.eq.1) then 164 clear=clear-springarea(2)*(1.d0-reltime) 165 endif 166 cstr(1)=clear 167! 168! 169 if(int(elcon(3,1,imat)).eq.1) then 170! 171! exponential overclosure 172! 173 if(dabs(elcon(2,1,imat)).lt.1.d-30) then 174 elas(1)=0.d0 175 beta=1.d0 176 else 177! 178 alpha=elcon(2,1,imat)*springarea(1) 179 beta=elcon(1,1,imat) 180 if(-beta*clear.gt.23.d0-dlog(alpha)) then 181 beta=(dlog(alpha)-23.d0)/clear 182 endif 183 elas(1)=dexp(-beta*clear+dlog(alpha)) 184 endif 185 if(nener.eq.1) then 186 senergy=elas(1)/beta; 187 endif 188 elseif((int(elcon(3,1,imat)).eq.2).or. 189 & (int(elcon(3,1,imat)).eq.4)) then 190! 191! linear overclosure/tied overclosure 192! 193 elas(1)=-springarea(1)*elcon(2,1,imat)*clear/kscale 194! 195! spring scaling for explicit dynamics 196! 197 if((mscalmethod.eq.2).or.(mscalmethod.eq.3)) then 198 elas(1)=elas(1)*smscale(nelem) 199 endif 200! 201 if(nener.eq.1) then 202 senergy=-elas(1)*clear/2.d0; 203 endif 204 elseif(int(elcon(3,1,imat)).eq.3) then 205! 206! tabular overclosure 207! 208! interpolating the material data 209! 210 call materialdata_sp(elcon,nelcon,imat,ntmat_,i,t1l, 211 & elconloc,kode,plicon,nplicon,npmat_,plconloc,ncmat_) 212 overlap=-clear 213 niso=int(plconloc(81)) 214 do i=1,niso 215 xiso(i)=plconloc(2*i-1) 216 yiso(i)=plconloc(2*i) 217 enddo 218 call ident(xiso,overlap,niso,id) 219 if(id.eq.0) then 220 pres=yiso(1) 221 if(nener.eq.1) then 222 senergy=yiso(1)*overlap; 223 endif 224 elseif(id.eq.niso) then 225 pres=yiso(niso) 226 if(nener.eq.1) then 227 senergy=yiso(1)*xiso(1) 228 do i=2,niso 229 senergy=senergy+(xiso(i)-xiso(i-1))* 230 & (yiso(i)+yiso(i-1))/2.d0 231 enddo 232 senergy=senergy+(pres-xiso(niso))*yiso(niso) 233 endif 234 else 235 xk=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id)) 236 pres=yiso(id)+xk*(overlap-xiso(id)) 237 if(nener.eq.1) then 238 senergy=yiso(1)*xiso(1) 239 do i=2, id 240 senergy=senergy+(xiso(i)-xiso(i-1))* 241 & (yiso(i)+yiso(i-1))/2.d0 242 enddo 243 senergy=senergy+(overlap-xiso(id))*(pres+yiso(id))/2.d0 244 endif 245 endif 246 elas(1)=springarea(1)*pres 247 if(nener.eq.1) senergy=springarea(1)*senergy 248 endif 249! 250! forces in the nodes of the contact element 251! 252 do i=1,3 253 fnl(i,nopep)=-elas(1)*xn(i) 254 enddo 255 cstr(4)=elas(1)/springarea(1) 256! 257! Coulomb friction for static calculations 258! 259 if((int(elcon(3,1,imat)).eq.4).and.(ncmat_.lt.7)) then 260 write(*,*) '*ERROR in springforc_f2f: for contact' 261 write(*,*) ' with pressure-overclosure=tied' 262 write(*,*) ' a stick slope using the *FRICTION' 263 write(*,*) ' card is mandatory' 264 call exit(201) 265 endif 266! 267 if(ncmat_.ge.7) then 268! 269! tied contact 270! 271 if(int(elcon(3,1,imat)).eq.4) then 272 um=1.d30 273 else 274 um=elcon(6,1,imat) 275 endif 276 stickslope=elcon(7,1,imat)/kscale 277! 278! spring scaling for explicit dynamics 279! 280 if((mscalmethod.eq.2).or.(mscalmethod.eq.3)) then 281 stickslope=stickslope*smscale(nelem) 282 endif 283! 284 if(um.gt.0.d0) then 285 if(1.d0 - dabs(xn(1)).lt.1.5231d-6) then 286! 287! calculating the local directions on master surface 288! 289 t1(1)=-xn(3)*xn(1) 290 t1(2)=-xn(3)*xn(2) 291 t1(3)=1.d0-xn(3)*xn(3) 292 else 293 t1(1)=1.d0-xn(1)*xn(1) 294 t1(2)=-xn(1)*xn(2) 295 t1(3)=-xn(1)*xn(3) 296 endif 297 dt1=dsqrt(t1(1)*t1(1)+t1(2)*t1(2)+t1(3)*t1(3)) 298 do i=1,3 299 t1(i)=t1(i)/dt1 300 enddo 301 t2(1)=xn(2)*t1(3)-xn(3)*t1(2) 302 t2(2)=xn(3)*t1(1)-xn(1)*t1(3) 303 t2(3)=xn(1)*t1(2)-xn(2)*t1(1) 304! 305! linear stiffness of the shear stress versus 306! slip curve 307! 308 xk=stickslope*springarea(1) 309! 310! calculating the relative displacement between the slave node 311! and its projection on the master surface 312! 313 do i=1,3 314 alnew(i)=vl(i,nopep) 315 do j=1,nopem 316 alnew(i)=alnew(i)-shp2m(4,j)*vl(i,j) 317 enddo 318 enddo 319! 320! calculating the difference in relative displacement since 321! the start of the increment = lamda^* 322! 323 do i=1,3 324 al(i)=alnew(i)-xstateini(3+i,1,ne0+igauss) 325 enddo 326! 327! ||lambda^*|| 328! 329 val=al(1)*xn(1)+al(2)*xn(2)+al(3)*xn(3) 330! 331! update the relative tangential displacement 332! 333 do i=1,3 334 t(i)=xstateini(6+i,1,ne0+igauss)+al(i)-val*xn(i) 335 enddo 336! 337! store the actual relative displacement and 338! the actual relative tangential displacement 339! 340 do i=1,3 341 xstate(3+i,1,ne0+igauss)=alnew(i) 342 xstate(6+i,1,ne0+igauss)=t(i) 343 enddo 344! 345! size of normal force 346! 347 dfnl=dsqrt(fnl(1,nopep)**2+fnl(2,nopep)**2+ 348 & fnl(3,nopep)**2) 349! 350! maximum size of shear force 351! 352 if(int(elcon(3,1,imat)).eq.4) then 353 dfshear=1.d30 354 else 355 dfshear=um*dfnl 356 endif 357! 358! plastic and elastic slip 359! 360 do i=1,3 361 tp(i)=xstateini(i,1,ne0+igauss) 362 te(i)=t(i)-tp(i) 363 enddo 364 dte=dsqrt(te(1)*te(1)+te(2)*te(2)+te(3)*te(3)) 365! 366! trial force 367! 368 do i=1,3 369 ftrial(i)=xk*te(i) 370 enddo 371 dftrial=xk*dte 372c dftrial=dsqrt(ftrial(1)**2+ftrial(2)**2+ftrial(3)**2) 373! 374! check whether stick or slip 375! 376 if((dftrial.lt.dfshear).or.(dftrial.le.0.d0).or. 377 & (ielas.eq.1)) then 378! 379! stick 380! 381 do i=1,3 382 fnl(i,nopep)=fnl(i,nopep)+ftrial(i) 383 xstate(i,1,ne0+igauss)=tp(i) 384 enddo 385 cstr(5)=(ftrial(1)*t1(1)+ftrial(2)*t1(2)+ 386 & ftrial(3)*t1(3))/springarea(1) 387 cstr(6)=(ftrial(1)*t2(1)+ftrial(2)*t2(2)+ 388 & ftrial(3)*t2(3))/springarea(1) 389! 390! shear elastic energy 391! 392 if(nener.eq.1) senergy=senergy+dftrial*dte 393 else 394! 395! slip 396! 397 dg=(dftrial-dfshear)/xk 398 do i=1,3 399 ftrial(i)=te(i)/dte 400 fnl(i,nopep)=fnl(i,nopep)+dfshear*ftrial(i) 401 xstate(i,1,ne0+igauss)=tp(i)+dg*ftrial(i) 402 enddo 403 cstr(5)=(dfshear*ftrial(1)*t1(1)+ 404 & dfshear*ftrial(2)*t1(2)+ 405 & dfshear*ftrial(3)*t1(3))/springarea(1) 406 cstr(6)=(dfshear*ftrial(1)*t2(1)+ 407 & dfshear*ftrial(2)*t2(2)+ 408 & dfshear*ftrial(3)*t2(3))/springarea(1) 409! 410! shear elastic and viscous energy 411! 412 if(nener.eq.1) then 413 senergy=senergy+dfshear*dfshear/xk 414 venergy=dg*dfshear 415 endif 416! 417 endif 418 endif 419! 420! storing the tangential displacements 421! 422 cstr(2)=t(1)*t1(1)+t(2)*t1(2)+t(3)*t1(3) 423 cstr(3)=t(1)*t2(1)+t(2)*t2(2)+t(3)*t2(3) 424 endif 425! 426! force in the master nodes 427! 428 do i=1,3 429 do j=1,nopem 430 fnl(i,j)=-shp2m(4,j)*fnl(i,nopep) 431 enddo 432 enddo 433! 434! forces in the nodes of the slave face 435! 436 do i=1,3 437 do j=1,nopes 438 fnl(i,nopem+j)=shp2s(4,j)*fnl(i,nopep) 439 enddo 440 enddo 441! 442! write statements for Malte Krack 443! 444c if(iout.gt.0) then 445c write(*,*) 'contact element: ',nelem 446c write(*,*) 'undeformed location of the integration point' 447c write(*,*) ((pl(j,nopep)-vl(j,nopep)),j=1,3) 448c write(*,*) 'deformed location of the integration point' 449c write(*,*) (pl(j,nopep),j=1,3) 450c write(*,*) 'nodes and shape values' 451c do j=1,nopem 452c write(*,*) konl(j),-shp2m(4,j) 453c enddo 454c do j=1,nopes 455c write(*,*) konl(nopem+j),shp2s(4,j) 456c enddo 457c write(*,*) 'contact force' 458c write(*,*) fnl(1,nopep),fnl(2,nopep),fnl(3,nopep) 459c write(*,*) 'slave area' 460c write(*,*) springarea(1) 461c endif 462! 463 return 464 end 465 466