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 checkimpacts(ne,neini,temax,sizemaxinc, 20 & energyref,tmin,tper,idivergence, 21 & iforceincsize,istab,dtheta,r_abs,energy,energyini,allwk, 22 & allwkini,dampwk,dampwkini,emax,mortar,maxdecay,enetoll) 23! 24! # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 25! Routine that contains the implementation of the logic to 26! rule the increment size during contact conditions. 27! The values of tolerances have been tested for the ball 28! model, the two sliding blades (simplified model of 29! blade from Mr. Wittig) and the real blade model. 30! Friction has not been tested deeply. 31! 32! Main variables and meaning 33! 34! sizemaxinc : maximum size of the current increment 35! iforceincsize : flag to set "dtheta=sizemaxinc" in 36! checkconvergence. 37! cstate : vector containing contact data 38! temax : max. natural period of oscillation 39! icase : flag to print debug informations 40! emax : maximum energy of the system over the ti- 41! me history 42! r : energy residual before contact (or re- 43! adapted) 44! delta : eneres normalized 45! fact : factor to set sizemaxinc according to the 46! contact formulation 47! stab_th : \hat{r}_{e}(t_n) -> mod belytschko before 48! contact (or initial value). This is the 49! stability threshold 50! delta_r_rel : \hat{r}_{e}(t) - \hat{r}_{e}(t-1) -> varia- 51! tion of the modified belitschko criterion 52! used to control jumps 53! r_rel : \hat{r}_{e}(t) actual value of the modified 54! belitschko criterion 55! 56! Proposed by Matteo Pacher 57! 58! # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 59! 60 implicit none 61! 62 integer idivergence, 63 & iforceincsize,ne,neini,istab,mortar 64! 65 real*8 temax,energyref,sizemaxinc,tmin,tper,dtheta, 66 & delta_r_rel,r_rel,delta,allwk,allwkini,energy(4), 67 & energyini(4),dampwk,dampwkini,emax,fact, 68 & r_rel_bc,maxdecay,r_abs,enetoll,delta_r_abs 69! 70! 71! 72! Initialization 73! 74 iforceincsize=0 75! 76! Adaption of the energy residual (absolute/relative check) 77! 78 if(dabs(r_abs).lt.(enetoll/4.d0))then 79 delta=r_abs*emax 80 else 81 delta=r_abs 82 endif 83! 84 if(mortar.eq.0)then 85 fact=10.d0 86 elseif(mortar.eq.1) then 87 fact=1.d0 88 endif 89! 90! Compute thresholds and energy values 91! 92 delta_r_abs=energy(1)+energy(2)+energy(3)+energy(4)-allwk- 93 & dampwk-(energyini(1)+energyini(2)+energyini(3)+ 94 & energyini(4)-allwkini-dampwkini) 95! 96 if(emax.le.0.d0)then 97! 98! No kinetic energy in the structure: energyref is the internal energy 99! this happens at the beginning of the calculation 100! 101 r_rel=(energy(1)+energy(2)+energy(3)+energy(4)-allwk- 102 & dampwk-energyref)/energyref 103 delta_r_rel=delta_r_abs/energyref 104 r_rel_bc=delta/energyref 105 else 106 r_rel=(energy(1)+energy(2)+energy(3)+energy(4)-allwk- 107 & dampwk-energyref)/emax 108 delta_r_rel=delta_r_abs/emax 109 r_rel_bc=delta/emax 110 endif 111! 112! Logic to adapt the increment size 113! 114 if(mortar.eq.0)then 115! 116! Energy conservation rules for NODE TO SURFACE penalty contact 117! 118 if((delta_r_rel.lt.(-0.008d0)).and.(ne.ge.neini))then 119! 120! Impact (or too high variation during pers. contact) 121! delta_r_rel = r_rel-r_rel_ini 122! 123 idivergence=1 124 sizemaxinc=dtheta*0.25d0 125 iforceincsize=1 126 elseif((r_rel-r_rel_bc.gt.0.0025d0).and.(ne.le.neini))then 127! 128! Rebound (or too high variation during pers. contact) 129! r_rel_bc is r_rel before contact 130! 131 idivergence=1 132 sizemaxinc=dtheta*0.5d0 133 iforceincsize=1 134 else 135! 136! Persistent Contact 137! 138 if(r_rel.gt.(-0.9d0*maxdecay))then 139 sizemaxinc=max(fact*temax/tper,1.01d0*dtheta) 140 sizemaxinc=min(sizemaxinc,100.d0*temax/tper) 141 else 142 sizemaxinc=max(temax/tper/10.d0,0.5d0*dtheta) 143 istab=1 144 endif 145! 146 endif 147! 148 elseif(mortar.eq.1)then 149! 150! Energy conservation rules for SURFACE TO SURFACE penalty contact 151! 152 if((delta_r_rel.lt.(-0.008d0)).and.(ne.ge.neini))then 153! 154! Impact (or too high variation during pers. contact) 155! delta_r_rel = r_rel-r_rel_ini 156! 157 idivergence=1 158 sizemaxinc=dtheta*0.25d0 159 iforceincsize=1 160! 161 elseif((r_rel-r_rel_bc.gt.0.0025d0).and.(ne.le.neini))then 162! 163! Rebound (or too high variation during pers. contact) 164! r_rel_bc is r_rel before contact 165! 166 idivergence=1 167 sizemaxinc=dtheta*0.5d0 168 iforceincsize=1 169! 170 else 171! 172! Persistent Contact 173! 174 if(r_rel.gt.(-0.9d0*maxdecay))then 175 sizemaxinc=min(fact*temax/tper,1.1d0*dtheta) 176 sizemaxinc=min(sizemaxinc,100.d0*temax/tper) 177 else 178 sizemaxinc=max(temax/tper/10.d0,0.5d0*dtheta) 179 istab=1 180 endif 181 endif 182 endif !(mortar) 183! 184 if(sizemaxinc.lt.tmin)then 185 sizemaxinc=tmin 186 endif 187! 188 write(*,*) '*INFO in checkimpacts: due to impact rules the' 189 write(*,*) ' maximum allowed time increment has been' 190 write(*,*) ' changed to',sizemaxinc*tper 191! 192 return 193 end 194