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 calcstabletimeinccont(ne,lakon,kon,ipkon,mi, 20 & ielmat,elcon,mortar,adb,alpha,nactdof,springarea, 21 & ne0,ntmat_,ncmat_,dtcont,smscale,dtset,mscalmethod) 22! 23! Calculates the critical time increment (CTI) for contact 24! spring elements based on the Courant 25! Criterion for Explicit Dynamics calculations. 26! 27! mscalmethod<0: no explicit dynamics 28! 0: explicit dynamics, no scaling 29! 1: explicit dynamics, volumetric scaling active 30! 2: explicit dynamics, contact scaling active 31! 3: explicit dynamics, volumetric and contact scaling 32! 33 implicit none 34! 35 character*8 lakon(*),lakonl 36! 37 integer j,ne,nope,kon(*),ipkon(*),indexe,nelem, 38 & ncmat_,ntmat_,mi(*),ielmat(mi(3),*),imat,mortar, 39 & indexn,nactdof(0:mi(2),*),ne0,nopem,mscalmethod, 40 & icount 41! 42 real*8 elcon(0:ncmat_,ntmat_,*),safefac,xk,adb(*),alpha(*), 43 & critom,damping,springms,springmm,springfac,dtcont,xmacont, 44 & springarea(2,*),areaslav,bet,gam,smscale(*),dtset 45! 46 dtcont=1.d30 47 safefac=0.80d0 48! 49 xmacont=0.0d0 50! 51 damping=0.d0 52 icount=0 53! 54 bet=(1.d0-alpha(1))*(1.d0-alpha(1))/4.d0 55 gam=0.5d0-alpha(1) 56! 57! Omega Critical 58! Om_cr=dt*freq_max 59! 60 critom=dsqrt(damping*damping*(1.d0+2.d0*alpha(1)*(1.d0-gam)) 61 & *(1.d0+2.d0*alpha(1)*(1.d0-gam)) 62 & + 2.d0*(gam+2.d0*alpha(1)*(gam-bet))) 63 critom=0.98d0*(-damping*(1.d0+2.d0*alpha(1)*(1.d0-gam))+critom) 64 & /(gam+2.d0*alpha(1)*(gam-bet)) !eq 25 miranda 65! 66! ** DO per element 67! 68 do nelem=ne0+1,ne 69 indexe=ipkon(nelem) 70 lakonl=lakon(nelem) 71 imat=ielmat(1,nelem) 72! 73 xk=elcon(2,1,imat) 74! 75 springmm=0.0d0 76 springms=0.0d0 77! 78 if(mortar.eq.0) then 79! 80! node-to-face 81! 82! nope is the total number of nodes: 83! master nodes+1 slave node 84! (notice -47 instead of -48) 85! 86 nope=ichar(lakonl(8:8))-47 87 springfac=1.d0 88! 89 do indexn=1,nope 90 xmacont=0.0d0 91! 92 do j=1,3 93 if(nactdof(j,kon(indexe+indexn)).gt.0)then 94 xmacont=max(xmacont, 95 & adb(nactdof(j,kon(indexe+indexn)))) 96 endif 97 enddo 98! 99 if(indexn.eq.nope)then 100 springms=springms+xmacont 101! 102! mean mass at master nodes 103! 104 springmm=springmm/(nope-1.0d0) 105 else 106 springmm=springmm+xmacont 107 endif 108 enddo 109! 110 areaslav=springarea(1,kon(indexe+nope+1)) 111 xk=xk*areaslav 112! 113! checking whether mass or slave side is constrained 114! 115 if((springmm.le.0.d0).and.(springms.le.0.d0)) then 116 cycle 117 elseif(springmm.le.0.d0) then 118 springmm=springms 119 elseif(springms.le.0.d0) then 120 springms=springmm 121 endif 122! 123 elseif(mortar.eq.1) then 124! 125! face-to-face 126! 127 nopem=ichar(lakonl(8:8))-48 128 springfac=0.1d0 129! 130 do indexn=1,kon(indexe) 131! 132 xmacont=0.0d0 133 do j=1,3 134 if(nactdof(j,kon(indexe+indexn)).gt.0)then 135 xmacont=max(xmacont, 136 & adb(nactdof(j,kon(indexe+indexn)))) 137 endif 138 enddo 139 if(indexn.gt.nopem)then !slave 140 springms=springms+xmacont 141 else !master 142 springmm=springmm+xmacont 143 endif 144 enddo 145! 146! mean mass at slave and master nodes 147! 148 springms=springms/(kon(indexe)-nopem) 149 springmm=springmm/nopem 150! 151 areaslav=springarea(1,kon(1+indexe+kon(indexe))) 152 xk=xk*areaslav 153! 154! checking whether mass or slave side is constrained 155! 156 if((springmm.le.0.d0).and.(springms.le.0.d0)) then 157 cycle 158 elseif(springmm.le.0.d0) then 159 springmm=springms 160 elseif(springms.le.0.d0) then 161 springms=springmm 162 endif 163 endif 164! 165! linear pressure-overclosure 166! 167 if(int(elcon(3,1,imat)).eq.2) then 168! 169 springmm=springmm/2.0d0 170 springms=springms/2.0d0 171! 172! time increment needed by contact spring 173! 174 smscale(nelem)=springfac*critom* 175 & dsqrt(springmm*springms/ 176 & ((springmm+springms)*xk)) 177! 178 dtcont = 179 & min(dtcont,smscale(nelem)) 180! 181 else 182 write(*,*) '*ERROR in calcstabletimeinccont:' 183 write(*,*) ' in explicit dynamic calculations' 184 write(*,*) ' only linear pressure-overclosure' 185 write(*,*) ' is allowed' 186 call exit(201) 187 endif 188 enddo 189! ** ENDDO per element 190! 191! check whether spring scaling is necessary in order to meet the 192! smallest allowed time increment specified by the user 193! 194 if(dtcont.lt.dtset/safefac) then 195 dtset=dtset/safefac 196 if(mscalmethod.eq.1) then 197 mscalmethod=3 198 else 199 mscalmethod=2 200 endif 201! 202 do nelem=ne0+1,ne 203 if(smscale(nelem).ge.dtset) then 204 smscale(nelem)=1.d0 205 else 206 smscale(nelem)=(smscale(nelem)/dtset)**2 207 icount=icount+1 208 endif 209 enddo 210! 211 write(*,*) 'Selective Spring Scaling is active' 212 write(*,*) 213 write(*,*) 'Scaling factor of time increment:',dtset/dtcont 214 write(*,*) 'Reduction of spring stiffness by (maximum) =', 215 & (dtset/dtcont)**2 216 write(*,*) 'In total ',icount,'elements were scaled of' 217 write(*,*) ' ',ne-ne0,' contact elements' 218 write(*,*) 219! 220 dtset=dtset*safefac 221! 222 else 223 if(mscalmethod.eq.3) then 224 mscalmethod=1 225 elseif(mscalmethod.eq.2) then 226 mscalmethod=0 227 endif 228 endif 229! 230 dtcont=dtcont*safefac 231! 232 return 233 end 234