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 mafillkcomp(nef,ipnei,neifa,neiel,vfa,xxn,area, 20 & au,ad,jq,irow,nzs,b,vel,umfa,alet,ale,gradkfa,xxi, 21 & body,volume,ielfa,lakonf,ifabou,nbody,neq, 22 & dtimef,velo,veloo,cvfa,hcfa,cvel,gradvel,xload,xrlfa, 23 & xxj,nactdohinv,a1,a2,a3,flux,nefa,nefb,iau6,xxni,xxnj, 24 & iturbulent,f1,of2,yy,umel,gradkel,gradoel,inlet,sc) 25! 26! filling the matrix for the conservation of energy 27! 28 implicit none 29! 30 logical knownflux 31! 32 character*8 lakonf(*) 33! 34 integer i,nef,indexf,ipnei(*),j,ifa,iel,neifa(*), 35 & neiel(*),jq(*),irow(*),nzs,ielfa(4,*),nefa,nefb, 36 & ipointer,ifabou(*),nbody,neq,indexb,nactdohinv(*), 37 & iau6(6,*),iturbulent,inlet(*) 38! 39 real*8 xflux,vfa(0:7,*),xxn(3,*),area(*),au(*),ad(*),b(neq), 40 & vel(nef,0:7),umfa(*),alet(*),ale(*),coef,gradkfa(3,*), 41 & xxi(3,*),body(0:3,*),volume(*),dtimef,velo(nef,0:7), 42 & veloo(nef,0:7),rhovol,cvel(*),gradvel(3,3,*),sk, 43 & cvfa(*),hcfa(*),div,xload(2,*),xrlfa(3,*),unt, 44 & xxj(3,*),a1,a2,a3,flux(*),xxnj(3,*),xxni(3,*),difcoef, 45 & f1(*),of2(*),yy(*),umel(*),gradkel(3,*),gradoel(3,*), 46 & cd,arg1,sc(*),constant 47! 48! 49! 50 do i=nefa,nefb 51! 52 if(iturbulent.eq.1) then 53! 54! k-epsilon model 55! 56 sk=1.d0 57 elseif(iturbulent.eq.2) then 58! 59! k-omega model 60! 61 sk=0.5d0 62 else 63! 64! BSL and SST model 65! 66 cd=max(2.d0*vel(i,5)*0.856d0* 67 & (gradkel(1,i)*gradoel(1,i)+ 68 & gradkel(2,i)*gradoel(2,i)+ 69 & gradkel(3,i)*gradoel(3,i))/vel(i,7), 70 & 1.d-20) 71 arg1=min(max(dsqrt(vel(i,6))/(0.09d0*vel(i,7)*yy(i)), 72 & 500.d0*umel(i)/(vel(i,5)*yy(i)**2*vel(i,7))), 73 & 4.d0*vel(i,5)*0.856d0*vel(i,6)/(cd*yy(i)**2)) 74 f1(i)=dtanh(arg1**4) 75 if(iturbulent.eq.3) then 76! 77! BSL model 78! 79 sk=0.5d0*f1(i)+(1.d0-f1(i)) 80 else 81! 82! SST model 83! 84 sk=0.85d0*f1(i)+(1.d0-f1(i)) 85 endif 86 endif 87! 88 do indexf=ipnei(i)+1,ipnei(i+1) 89! 90! convection 91! 92 ifa=neifa(indexf) 93 iel=neiel(indexf) 94 xflux=flux(indexf) 95! 96 if(xflux.ge.0.d0) then 97! 98! outflowing flux 99! 100 ad(i)=ad(i)+xflux 101! 102 b(i)=b(i)-(vfa(6,ifa)-vel(i,6))*xflux 103! 104 else 105 if(iel.gt.0) then 106! 107! incoming flux from neighboring element 108! 109 au(indexf)=au(indexf)+xflux 110! 111 b(i)=b(i)-(vfa(6,ifa)-vel(iel,6))*xflux 112! 113 else 114! 115! incoming flux through boundary 116! 117 if(ielfa(2,ifa).lt.0) then 118 if(inlet(ifa).eq.1) then 119 b(i)=b(i)-vfa(6,ifa)*xflux 120 endif 121 endif 122 endif 123 endif 124! 125! diffusion 126! 127 if(iturbulent.le.3) then 128! 129! k-epsilon, k-omega or BSL model 130! 131 difcoef=umfa(ifa)+sk*vfa(5,ifa)*vfa(6,ifa)/vfa(7,ifa) 132 else 133! 134! SST model 135! 136 difcoef=umfa(ifa)+sk*vfa(5,ifa)*(0.31d0*vfa(6,ifa))/ 137 & max(0.31d0*vfa(7,ifa),of2(i)) 138 endif 139! 140 if(iel.ne.0) then 141! 142! neighboring element 143! 144 coef=difcoef*alet(indexf) 145 ad(i)=ad(i)+coef 146 au(indexf)=au(indexf)-coef 147! 148! correction for non-orthogonal grid 149! 150 b(i)=b(i)+difcoef* 151 & (gradkfa(1,ifa)*xxnj(1,indexf)+ 152 & gradkfa(2,ifa)*xxnj(2,indexf)+ 153 & gradkfa(3,ifa)*xxnj(3,indexf)) 154 else 155! 156! boundary 157! 158 knownflux=.false. 159 ipointer=abs(ielfa(2,ifa)) 160 if(ipointer.gt.0) then 161 if(ielfa(3,ifa).gt.0) then 162! 163! if diffusion is not necessarily equal to zero: 164! turbulent variable is assumed fixed 165! 166 coef=difcoef*ale(indexf) 167 ad(i)=ad(i)+coef 168 b(i)=b(i)+coef*vfa(6,ifa) 169 else 170! 171! outlet: no diffusion 172! 173 endif 174 endif 175! 176! correction for non-orthogonal grid 177! 178 if(.not.knownflux) then 179 b(i)=b(i)+difcoef* 180 & (gradkfa(1,ifa)*xxnj(1,indexf)+ 181 & gradkfa(2,ifa)*xxnj(2,indexf)+ 182 & gradkfa(3,ifa)*xxnj(3,indexf)) 183 endif 184 endif 185 enddo 186! 187 rhovol=vel(i,5)*volume(i) 188 div=gradvel(1,1,i)+gradvel(2,2,i)+gradvel(3,3,i) 189! 190! sink terms are treated implicitly (lhs) 191! 192 ad(i)=ad(i)+rhovol*0.09d0*vel(i,7) 193! 194! production term containing the first invariant of the 195! Reynolds stress (= the turbulent kinetic energy) 196! 197 if(div.gt.0.d0) then 198 ad(i)=ad(i)+2.d0*rhovol*div/3.d0 199 else 200 b(i)=b(i)-2.d0*rhovol*vel(i,6)*div/3.d0 201 endif 202! 203! rest of the production term 204! (positive part is treated explicitly (rhs), 205! negative part is treated implicitly (lhs)) 206! 207 if(iturbulent.le.3) then 208! 209! k-epsilon, k-omega and BSL-model: the turbulent 210! viscosity=k/omega 211! 212 unt=vel(i,6)/vel(i,7) 213 b(i)=b(i)+rhovol*unt*( 214 & (2.d0*(gradvel(1,1,i)**2+gradvel(2,2,i)**2+ 215 & gradvel(3,3,i)**2)+ 216 & (gradvel(1,2,i)+gradvel(2,1,i))**2+ 217 & (gradvel(1,3,i)+gradvel(3,1,i))**2+ 218 & (gradvel(2,3,i)+gradvel(3,2,i))**2)) 219 ad(i)=ad(i)+2.d0*rhovol*div**2/(3.d0*vel(i,7)) 220 else 221! 222! SST model: other definition of the turbulent 223! viscosity 224! 225 unt=(0.31d0*vel(i,6))/max(0.31d0*vel(i,7),of2(i)) 226 b(i)=b(i)+rhovol*unt*( 227 & (2.d0*(gradvel(1,1,i)**2+gradvel(2,2,i)**2+ 228 & gradvel(3,3,i)**2)+ 229 & (gradvel(1,2,i)+gradvel(2,1,i))**2+ 230 & (gradvel(1,3,i)+gradvel(3,1,i))**2+ 231 & (gradvel(2,3,i)+gradvel(3,2,i))**2)) 232 ad(i)=ad(i)+2.d0*rhovol*div**2/(3.d0*vel(i,7)) 233 endif 234! 235! transient term 236! 237 constant=rhovol*sc(i) 238 b(i)=b(i)-(a2*velo(i,6)+a3*veloo(i,6))*constant 239 rhovol=a1*constant 240 ad(i)=ad(i)+constant 241! 242 enddo 243! 244 return 245 end 246