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 calcgammak(nface,ielfa,vel,gradkel,gamma,xlet, 20 & xxn,xxj,ipnei,betam,nef,flux) 21! 22! determine gamma for the turbulent kinetic energy: 23! upwind difference: gamma=0 24! central difference: gamma=1 25! 26 implicit none 27! 28 integer nface,ielfa(4,*),i,j,indexf,ipnei(*),iel1,iel2,nef 29! 30 real*8 vel(nef,0:7),gradkel(3,*),xxn(3,*),xxj(3,*),vud,vcd, 31 & gamma(*),phic,xlet(*),betam,flux(*) 32! 33 do i=1,nface 34 iel2=ielfa(2,i) 35! 36! faces with only one neighbor need not be treated 37! 38 if(iel2.le.0) cycle 39 iel1=ielfa(1,i) 40 j=ielfa(4,i) 41 indexf=ipnei(iel1)+j 42! 43 vcd=vel(iel2,6)-vel(iel1,6) 44 if(dabs(vcd).lt.1.d-3*dabs(vel(iel1,6))) vcd=0.d0 45! 46 if(flux(indexf).ge.0.d0) then 47 vud=2.d0*xlet(indexf)* 48 & (gradkel(1,iel1)*xxj(1,indexf)+ 49 & gradkel(2,iel1)*xxj(2,indexf)+ 50 & gradkel(3,iel1)*xxj(3,indexf)) 51 else 52 vud=2.d0*xlet(indexf)* 53 & (gradkel(1,iel2)*xxj(1,indexf)+ 54 & gradkel(2,iel2)*xxj(2,indexf)+ 55 & gradkel(3,iel2)*xxj(3,indexf)) 56 endif 57! 58 if(dabs(vud).lt.1.d-20) then 59 gamma(i)=0.d0 60 cycle 61 endif 62! 63 phic=1.d0-vcd/vud 64! 65 if(phic.ge.1.d0) then 66 gamma(i)=0.d0 67 elseif(phic.le.0.d0) then 68 gamma(i)=0.d0 69 elseif(betam.le.phic) then 70 gamma(i)=1.d0 71 else 72 gamma(i)=phic/betam 73 endif 74c write(*,*) 'calcgammat',i,gamma(i) 75 enddo 76! 77 return 78 end 79