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 calcstabletimeincvol(ne0,elcon,nelcon, 20 & rhcon,nrhcon,alcon,nalcon,orab,ntmat_,ithermal,alzero, 21 & plicon,nplicon,plkcon,nplkcon,npmat_,mi,dtime, 22 & xstiff,ncmat_,vold,ielmat,t0,t1, 23 & matname,lakon,wavespeed,nmat,ipkon,co,kon,dtvol,alpha, 24 & smscale,dtset,mscalmethod) 25! 26! ************** 27! ------------Wavespeed calculation------------------------CARLO MT 28! Calculates the propagation wave speed in a material, selecting 29! appropiate procedure between isotropic, single crystals, and 30! anisotropic materials. All other cases of orthotropy are treated 31! as anisotropic. 32 33! Based on the procedure in: 34! C. Lane. The Development of a 2D Ultrasonic Array Inspection 35! for Single Crystal Turbine Blades. 36! Switzerland: Springer International Publishing, 2014. 37 38! ------------Critical time step calculation---------------CARLO MT 39! Calculates the critical time increment (CTI) based on the Courant 40! Criterion for Explicit Dynamics calculations. Temperature is 41! assumed averaged from the centroid of the element and material 42! wave propagation speeds must be calculated before. 43 44! ------------Mass Scaling -----------------------------CATHARINA C 45! 46! 47! mscalmethod<0: no explicit dynamics 48! 0: explicit dynamics, no scaling 49! 1: explicit dynamics, volumetric scaling active 50! 2: explicit dynamics, contact scaling active 51! 3: explicit dynamics, volumetric and contact scaling 52! 53 implicit none 54! 55 character*80 matname(*),amat 56 character*8 lakon(*),lakonl 57! 58 integer i,j,i1,nelem,ne0,nelcon(2,*),nrhcon(*),nalcon(2,*),imat, 59 & ntmat_,ithermal(*),mattyp,ihyper,istiff,kode,mi(*),kk, 60 & nplicon(0:ntmat_,*),nplkcon(0:ntmat_,*),ncmat_,iorth, 61 & ielmat(mi(3),*),nope,iorien,ipkon(*),null, 62 & konl(26),nopered,npmat_,nmat,kon(*),indexe,iflag,nopes, 63 & nfaces,ig,ifaceq(8,6),ifacet(6,4),ifacew(8,5),ielemmin, 64 & mscalmethod,icount 65! 66 real*8 elas(21),wavespeed(*),rhcon(0:1,ntmat_,*) , 67 & alcon(0:6,ntmat_,*),coords(3),orab(7,*),rho,alzero(*), 68 & t0l,t1l,elconloc(21),eth(6),plicon(0:2*npmat_,ntmat_,*), 69 & plkcon(0:2*npmat_,ntmat_,*),plconloc(802),dtime, 70 & xstiff(27,mi(1),*),elcon(0:ncmat_,ntmat_,*), 71 & t0(*),t1(*),shp(4,26),vold(0:mi(2),*),tt, 72 & e,un,um,al,wavspd,xi,et,ze,weight,co(3,*),xl(3,26),xsj, 73 & xl2(3,9),xsj2(3),xs2(3,7),shp2(7,9),hmin,area, 74 & volume,dtvol,safefac,alpha,bet,gam,critom,damping, 75 & geomfac,quadfac,smscale(*),dtset 76! 77 data ifaceq /4,3,2,1,11,10,9,12, 78 & 5,6,7,8,13,14,15,16, 79 & 1,2,6,5,9,18,13,17, 80 & 2,3,7,6,10,19,14,18, 81 & 3,4,8,7,11,20,15,19, 82 & 4,1,5,8,12,17,16,20/ 83 data ifacet /1,3,2,7,6,5, 84 & 1,2,4,5,9,8, 85 & 2,3,4,6,10,9, 86 & 1,4,3,8,10,7/ 87 data ifacew /1,3,2,9,8,7,0,0, 88 & 4,5,6,10,11,12,0,0, 89 & 1,2,5,4,7,14,10,13, 90 & 2,3,6,5,8,15,11,14, 91 & 4,6,3,1,12,15,9,13/ 92! 93 include "gauss.f" 94! 95 iflag=2 96 dtvol=1.d30 97 safefac=0.80d0 98 quadfac=0.3d0 99! 100 damping=0 101 icount=0 102! 103 bet=(1.d0-alpha)*(1.d0-alpha)/4.d0 104 gam=0.5d0-alpha 105! 106! Calcualtion of Omega Critical 107! Om_cr=dt*freq_max 108 critom=dsqrt(damping*damping*(1.d0+2.d0*alpha*(1.d0-gam)) 109 & *(1.d0+2.d0*alpha*(1.d0-gam)) 110 & +2.d0*(gam+2.d0*alpha*(gam-bet)) ) 111 critom=0.98d0*(-damping*(1.d0+2.d0*alpha*(1.d0-gam))+critom) 112 & /(gam+2.d0*alpha*(gam-bet)) !eq 25 miranda 113 114 write(*,*)'++CMT: Calculating Material Wave Speeds...' 115 write(*,*) 116! 117 null=0 118! 119! Initialization of wavespeed 120! 121 do i=1,nmat 122 wavespeed(i)=-1.d0 123 enddo 124! 125! ** DO per element 126 do nelem=1,ne0 127 if(ipkon(nelem).lt.0) cycle 128! 129 lakonl=lakon(nelem) 130 imat=ielmat(1,nelem) 131 amat=matname(imat) 132 wavspd=wavespeed(imat) 133 indexe=ipkon(nelem) 134! 135 if(lakonl(1:2).ne.'ES') then 136! 137! orientation is not important for the wave speed 138 iorien=0 139! 140 if(nelcon(1,imat).lt.0) then 141 ihyper=1 142 else 143 ihyper=0 144 endif 145! 146! ------------Shape Functions -------------------------------- 147! Element type , Missing: C3D10T, C3D8R, C3D20R ? 148 ! 149 geomfac=1.d0 150 ! 151 if(lakon(nelem)(4:5).eq.'20')then 152 nope=20 153 nopes=8 154 nfaces=6 155 geomfac=quadfac 156 elseif(lakon(nelem)(1:5).eq.'C3D8I')then 157 nope=8 158 nopes=4 159 nfaces=6 160 geomfac=quadfac 161 geomfac=0.5d0 162 elseif(lakon(nelem)(4:4).eq.'8') then 163 nope=8 164 nopes=4 165 nfaces=6 166 elseif(lakon(nelem)(4:5).eq.'10')then 167 nope=10 168 nopes=6 169 nfaces=4 170 geomfac=quadfac 171 elseif(lakon(nelem)(4:4).eq.'4') then 172 nope=4 173 nopes=3 174 nfaces=4 175 elseif(lakon(nelem)(4:5).eq.'15')then 176 nope=15 177 nfaces=5 178 geomfac=quadfac 179 elseif(lakon(nelem)(4:4).eq.'6') then 180 nope=6 181 nfaces=5 182 else 183 cycle 184 endif 185! 186! Find center of the element for avg temp value on the element to 187! get properties later 188! if HEX 189 if((lakon(nelem)(4:5).eq.'20').or. 190 & (lakon(nelem)(4:4).eq.'8')) then 191 xi=0.d0 192 et=0.d0 193 ze=0.d0 194 weight=8.d0 195! if TET 196 elseif((lakon(nelem)(4:5).eq.'10').or. 197 & (lakon(nelem)(4:4).eq.'4')) then 198 xi=gauss3d4(1,1) !all integration points mint3d 199 et=gauss3d4(2,1) 200 ze=gauss3d4(3,1) 201 weight=weight3d4(1) 202! elseif WEDGES 203 elseif((lakonl(4:5).eq.'15').or. 204 & (lakonl(4:4).eq.'6'))then 205 xi=1.d0/3.d0 206 et=1.d0/3.d0 207 ze=0.d0 208 weight=1.d0 209 endif 210! 211! computation of the coordinates of the local nodes 212 do i=1,nope 213 konl(i)=kon(indexe+i) 214 do j=1,3 215 xl(j,i)=co(j,konl(i)) 216 enddo 217 enddo 218 ! 219 if (nope.eq.20)then 220 call shape20h(xi,et,ze,xl,xsj,shp,iflag) 221 elseif(nope.eq.8) then 222 call shape8h(xi,et,ze,xl,xsj,shp,iflag) 223 elseif(nope.eq.10)then 224 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 225 elseif(nope.eq.4) then 226 call shape4tet (xi,et,ze,xl,xsj,shp,iflag) 227 elseif(nope.eq.15)then 228 call shape15w(xi,et,ze,xl,xsj,shp,iflag) 229 else 230 call shape6w(xi,et,ze,xl,xsj,shp,iflag) 231 endif 232! 233! ------------Wavespeed calcualtion-------------------------------- 234! calculating the temperature in the integration point 235! 236 kk=1 ! Only 1 integration point is considered, CENTER 237 t0l=0.d0 238 t1l=0.d0 239 if(ithermal(1).eq.1) then 240 if(lakonl(4:5).eq.'8 ') then 241 do i1=1,nope 242 t0l=t0l+t0(konl(i1))/8.d0 243 t1l=t1l+t1(konl(i1))/8.d0 244 enddo 245 elseif(lakonl(4:6).eq.'20 ')then 246 nopered=20 247 call lintemp(t0,konl,nopered,kk,t0l) 248 call lintemp(t1,konl,nopered,kk,t1l) 249 elseif(lakonl(4:6).eq.'10T') then 250 call linscal10(t0,konl,t0l,null,shp) 251 call linscal10(t1,konl,t1l,null,shp) 252 else 253 do i1=1,nope 254 t0l=t0l+shp(4,i1)*t0(konl(i1)) 255 t1l=t1l+shp(4,i1)*t1(konl(i1)) 256 enddo 257 endif 258 elseif(ithermal(1).ge.2) then 259 if(lakonl(4:5).eq.'8 ') then 260 do i1=1,nope 261 t0l=t0l+t0(konl(i1))/8.d0 262 t1l=t1l+vold(0,konl(i1))/8.d0 263 enddo 264 elseif(lakonl(4:6).eq.'20 ')then 265 nopered=20 266 call lintemp_th0(t0,konl,nopered,kk,t0l,mi) 267 call lintemp_th1(vold,konl,nopered,kk,t1l,mi) 268 elseif(lakonl(4:6).eq.'10T') then 269 call linscal10(t0,konl,t0l,null,shp) 270 call linscal10(vold,konl,t1l,mi(2),shp) 271 else 272 do i1=1,nope 273 t0l=t0l+shp(4,i1)*t0(konl(i1)) 274 t1l=t1l+shp(4,i1)*vold(0,konl(i1)) 275 enddo 276 endif 277 endif 278 tt=t1l-t0l 279! 280 kode=nelcon(1,imat) 281! 282! material data 283 istiff=1 284 call materialdata_me(elcon,nelcon,rhcon,nrhcon,alcon,nalcon, 285 & imat,amat,iorien,coords,orab,ntmat_,elas,rho, 286 & nelem,ithermal,alzero,mattyp,t0l,t1l, 287 & ihyper,istiff,elconloc,eth,kode,plicon, 288 & nplicon,plkcon,nplkcon,npmat_, 289 & plconloc,mi(1),dtime,kk, 290 & xstiff,ncmat_) 291! 292 if(mattyp.eq.1) then 293 e=elas(1) 294 un=elas(2) 295 um=e/(1.d0+un) 296 al=un*um/(1.d0-2.d0*un) 297 um=um/2.d0 298 wavspd=max(wavspd, 299 & dsqrt((e*(1-un))/((1+un)*(1-2*un)*rho))) 300 elseif(mattyp.eq.2) then 301! 302! single crystal 303 if(((elas(1).eq.elas(3)).and.(elas(1).eq.elas(6)).and. 304 & (elas(3).eq.elas(6))).and. 305 & ((elas(2).eq.elas(4)).and.(elas(2).eq.elas(5)).and. 306 & (elas(4).eq.elas(5))).and. 307 & ((elas(7).eq.elas(8)).and.(elas(7).eq.elas(9)).and. 308 & (elas(8).eq.elas(9)))) then 309 wavspd=max(wavspd,dsqrt((1/3.d0)*(elas(1)+2.0*elas(2)+ 310 & 4.0d0*elas(7))/rho)) 311 wavspd=max(wavspd,dsqrt((1/2.d0)* 312 & (elas(1)+elas(2)+2.0*elas(7))/rho)) 313 wavspd=max(wavspd,dsqrt(elas(1)/rho)) 314 else 315 iorth=1 316 call anisomaxwavspd(elas,rho,iorth,wavspd ) 317 endif 318 elseif(mattyp.eq.3) then 319 iorth=0 320 call anisomaxwavspd(elas,rho,iorth,wavspd) 321 endif 322! 323 wavespeed(imat)=wavspd 324! 325! ------------Critical time step calcualtion----------------------- 326! 327! Divides volume accordingly per geometry of element 328! Carlo MT proposal 329! if HEX 330 if((lakon(nelem)(4:5).eq.'20').or. 331 & (lakon(nelem)(4:4).eq.'8')) then 332 volume=weight*xsj 333! if TET 334 elseif((lakon(nelem)(4:5).eq.'10').or. 335 & (lakon(nelem)(4:4).eq.'4')) then 336 volume=weight*xsj/3.d0 337! if WEDGES 338 elseif ( (lakonl(4:5).eq.'15').or. 339 & (lakonl(4:4).eq.'6'))then 340 volume=weight*xsj/2.d0 341 endif 342! 343 hmin=1.d30 344! 345! DO over sides 346 do ig=1,nfaces 347 if(lakon(nelem)(4:4).eq.'6')then 348 if(ig.le.2)then 349 nopes=3 350 else 351 nopes=4 352 endif 353 elseif(lakon(nelem)(4:5).eq.'15')then 354 if(ig.le.2)then 355 nopes=6 356 else 357 nopes=8 358 endif 359 endif 360 ! 361 if((nope.eq.20).or.(nope.eq.8))then 362 do i=1,nopes 363 do j=1,3 364 xl2(j,i)=co(j,konl(ifaceq(i,ig))) 365 enddo 366 enddo 367 elseif((nope.eq.10).or.(nope.eq.4))then 368 do i=1,nopes 369 do j=1,3 370 xl2(j,i)=co(j,konl(ifacet(i,ig))) 371 enddo 372 enddo 373 else 374 do i=1,nopes 375 do j=1,3 376 xl2(j,i)=co(j,konl(ifacew(i,ig))) 377 enddo 378 enddo 379 endif 380! 381 if((nopes.eq.4).or.(nopes.eq.8))then 382 xi=0.d0 383 et=0.d0 384 weight=4.d0 385 else 386 xi=1.d0/3.d0 387 et=1.d0/3.d0 388 weight=0.5d0 389 endif 390! 391 if (nopes.eq.8) then 392 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 393 elseif(nopes.eq.4) then 394 call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag) 395 elseif(nopes.eq.6) then 396 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 397c elseif(nopes.eq.7) then 398c call shape7tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 399 else 400 call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 401 endif 402! 403 area=weight*dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 404 & xsj2(3)*xsj2(3)) 405 406 hmin=min(hmin,(volume/area)) 407! 408 enddo 409! ENDDO over sides 410! 411! smallest element 412 if(critom/2*hmin/wavspd*geomfac.lt.dtvol)then 413 ielemmin=nelem 414 endif 415! 416! scaling of element: time increment required by element 417! 418 smscale(nelem)=critom/2*hmin/wavspd*geomfac 419! 420! smallest dtvol 421 dtvol=min(dtvol,smscale(nelem)) 422 endif !endif(lakonl(1:2).ne.'ES') 423! 424 enddo 425! ** ENDDO per element 426! 427 do i=1,nmat 428 write(*,*) 'Wave Speed for mat. ',i,wavespeed(i) 429 enddo 430 write(*,*) 431! 432! ------------Mass Scaling ---------------------------------------- 433! mscalmethod=1: selective mass scaling SMS 434! 435 if(dtvol.lt.dtset/safefac)then 436 dtset=dtset/safefac 437 mscalmethod=1 438! 439 open(40,file='WarnElementMassScaled.nam',status='unknown') 440 write(40,*) '*ELSET,ELSET=MassScaled' 441! 442 do nelem=1,ne0 443! scaling matrix 444 if(smscale(nelem).ge.dtset)then 445 smscale(nelem)=0.d0 446 else 447 smscale(nelem)=((dtset/smscale(nelem))**2)-1!beta 448 icount=icount+1 449 write(40,*) nelem 450 endif 451 enddo 452! 453 if(icount.gt.0) then 454 write(*,*) 'Selective Mass Scaling is active' 455 write(*,*) 456 write(*,*) 'Scaling factor of time increment:',dtset/dtvol 457 write(*,*) 'Overall mass is not changed:' 458 write(*,*) 'Manipulation of M matrix by beta (maximum) =', 459 & (((dtset/dtvol)**2)-1) 460 write(*,*) 'In total ',icount,'elements were scaled' 461 write(*,*) 462 write(*,*) '*INFO in calcstabletimeincvol:' 463 write(*,*) ' scaled nodes are stored in file' 464 write(*,*) ' WarnElementMassScaled.nam' 465 write(*,*) ' This file can be loaded into' 466 write(*,*) ' an active cgx-session by typing' 467 write(*,*) 468 & ' read WarnElementMassScaled.nam inp' 469 write(*,*) 470 close(40) 471 else 472 close(40,status='delete') 473 endif 474! 475 dtset=dtset*safefac 476 endif !endif (dtvol.lt.dtset) 477! 478 dtvol=dtvol*safefac 479! 480 return 481 end 482