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 compdt(nk,dt,nshcon,shcon,vold,ntmat_, 20 & iponoel,inoel,dtimef,ielmat,dh,cocon, 21 & ncocon,ithermal,mi, 22 & vcon,compressible,tincf,ierr,ifreesurface,dgravity,iit) 23! 24! - determine the time step for each node (stored in field dt 25! and the minimum value across all nodes (dtimef) 26! 27 implicit none 28! 29 integer nk,i,iponoel(*),inoel(2,*),index,nelem,ithermal(*), 30 & mi(*),compressible,ierr,ifreesurface,iit, 31 & nshcon(*),ntmat_,ielmat(mi(3),*),imat,ncocon(2,*), 32 & iflag 33! 34 real*8 dtimef,dt(*),dvi,r,cp,rho,shcon(0:3,ntmat_,*),tincf, 35 & vold(0:mi(2),*),temp,vel,dtcon,dtmed,dtimefold, 36 & dh(*),cocon(0:6,ntmat_,*),dtthd,cond,c, 37 & vcon(nk,0:mi(2)),dgravity 38! 39 data iflag /3/ 40! 41! determining the time increment dt for each node. 42! 43! edge nodes (fields iponoel and inoel are determined in precfd.f) 44! 45 dtimefold=dtimef 46 dtimef=1.d30 47! 48 do i=1,nk 49 index=iponoel(i) 50 if(index.le.0) cycle 51! 52! look at an element belonging to the edge node 53! 54 nelem=inoel(1,index) 55! 56! determining the time increment 57! 58 imat=ielmat(1,nelem) 59 temp=vold(0,i) 60! 61! convective time step (dtcon) 62! 63 vel=dsqrt(vold(1,i)**2+vold(2,i)**2+vold(3,i)**2) 64! 65 if(compressible.eq.1) then 66 if(ifreesurface.eq.0) then 67! 68! gas: speed of sound = dsqrt(kappa*r*temp) 69! 70 rho=vcon(i,4) 71 call materialdata_cp(imat,ntmat_,temp,shcon,nshcon,cp) 72 r=shcon(3,1,imat) 73 dtcon=dh(i)/(dsqrt(cp*r*temp/(cp-r))+vel) 74 c=cp*r*temp/(cp-r) 75 else 76! 77! shallow water: speed of sound = dsqrt(g*depth) 78! 79 dtcon=dh(i)/(dsqrt(dgravity*vcon(i,4))+vel) 80 endif 81 else 82! 83! liquid 84! 85 rho=vcon(i,4) 86 if(vel.lt.1.d-10) vel=1.d-10 87 dtcon=dh(i)/vel 88 endif 89 dt(i)=dtcon 90! 91! mechanical diffusion time step (dtmed) 92! 93 call materialdata_dvifem(imat,ntmat_,temp,shcon,nshcon,dvi) 94 if(dvi.gt.1.d-20) then 95 dtmed=dh(i)*dh(i)*rho/(2.d0*dvi) 96 dt(i)=dtcon*dtmed/(dtcon+dtmed) 97 endif 98! 99! thermal diffusion time step (dtthd) 100! 101 if(ithermal(1).gt.1) then 102 call materialdata_cond(imat,ntmat_,temp,cocon,ncocon, 103 & cond) 104 call materialdata_cp(imat,ntmat_,temp,shcon,nshcon,cp) 105 if(cond.gt.1.d-20) then 106 dtthd=dh(i)*dh(i)*rho*cp/(2.d0*cond) 107 dt(i)=(dt(i)*dtthd)/(dt(i)+dtthd) 108 endif 109 endif 110! 111! safety factor (default for compressible fluids: 1.25, 112! default for incompressible fluids: 1.00; 113! can be changed by user: fifth entry underneath *CFD) 114! 115 dt(i)=dt(i)/tincf 116! 117 if(dt(i).lt.dtimef) then 118 dtimef=dt(i) 119 endif 120! 121 enddo 122! 123! increased damping for incompressible fluids 124! the use of an internal (for the damping) and an external 125! (for the time derivative) time step stems from Zienkiewicz, 126! Taylor and Nithiarasu, The Finite Element Method for Fluid 127! Dynamics, 6th edition, p94 bottom and p 95 top. 128! 129 write(*,*) 130 write(*,*) 'dtimef= ',dtimef 131 write(*,*) 132 if(dtimef.le.0.d0) then 133 write(*,*) '*ERROR in compdt;' 134 write(*,*) ' negative time increment;' 135 write(*,*) ' the solution diverged' 136 ierr=1 137 elseif(((iit.gt.1).and.(dtimef.lt.dtimefold/100.d0))) then 138 write(*,*) '*ERROR in compdt;' 139 write(*,*) ' strongly decreasing time increment;' 140 write(*,*) ' the solution diverged' 141 ierr=1 142 endif 143! 144 return 145 end 146