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 mafillk(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,gamma,xrlfa, 23 & xxj,nactdohinv,a1,a2,a3,flux,nefa,nefb,iau6,xxni,xxnj, 24 & iturbulent,f1,of2,yy,umel,gradkel,gradoel,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 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,*),gamma(*),xrlfa(3,*), 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)-gamma(ifa)*(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)-gamma(ifa)*(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 indexb=-ielfa(2,ifa) 119 if(((ifabou(indexb+1).ne.0).and. 120 & (ifabou(indexb+2).ne.0).and. 121 & (ifabou(indexb+3).ne.0)).or. 122 & (dabs(xflux).lt.1.d-10)) then 123 b(i)=b(i)-vfa(6,ifa)*xflux 124 endif 125 endif 126 endif 127 endif 128! 129! diffusion 130! 131 if(iturbulent.le.3) then 132! 133! k-epsilon, k-omega or BSL model 134! 135 difcoef=umfa(ifa)+sk*vfa(5,ifa)*vfa(6,ifa)/vfa(7,ifa) 136 else 137! 138! SST model 139! 140 difcoef=umfa(ifa)+sk*vfa(5,ifa)*(0.31d0*vfa(6,ifa))/ 141 & max(0.31d0*vfa(7,ifa),of2(i)) 142 endif 143! 144 if(iel.ne.0) then 145! 146! neighboring element 147! 148 coef=difcoef*alet(indexf) 149 ad(i)=ad(i)+coef 150 au(indexf)=au(indexf)-coef 151! 152! correction for non-orthogonal grid 153! 154 b(i)=b(i)+difcoef* 155 & (gradkfa(1,ifa)*xxnj(1,indexf)+ 156 & gradkfa(2,ifa)*xxnj(2,indexf)+ 157 & gradkfa(3,ifa)*xxnj(3,indexf)) 158 else 159! 160! boundary; either temperature given or adiabatic 161! or outlet 162! 163 knownflux=.false. 164 ipointer=abs(ielfa(2,ifa)) 165 if(ipointer.gt.0) then 166 if((ifabou(ipointer+5).ne.0).or. 167 & (ifabou(ipointer+1).ne.0).or. 168 & (ifabou(ipointer+2).ne.0).or. 169 & (ifabou(ipointer+3).ne.0)) then 170! 171! no outlet: 172! (i.e. no wall || no sliding || at least one velocity given) 173! turbulent variable is assumed fixed 174! 175 coef=difcoef*ale(indexf) 176 ad(i)=ad(i)+coef 177 b(i)=b(i)+coef*vfa(6,ifa) 178 else 179! 180! outlet: no diffusion 181! 182 endif 183 endif 184! 185! correction for non-orthogonal grid 186! 187 if(.not.knownflux) then 188 b(i)=b(i)+difcoef* 189 & (gradkfa(1,ifa)*xxnj(1,indexf)+ 190 & gradkfa(2,ifa)*xxnj(2,indexf)+ 191 & gradkfa(3,ifa)*xxnj(3,indexf)) 192 endif 193 endif 194 enddo 195! 196! viscous dissipation 197! 198 rhovol=vel(i,5)*volume(i) 199! 200! sink terms are treated implicitly (lhs) 201! 202 ad(i)=ad(i)+rhovol*0.09d0*vel(i,7) 203! 204! source terms are treated explicitly (rhs) 205! 206 if(iturbulent.le.3) then 207! 208! k-epsilon, k-omega and BSL-model: the turbulent 209! viscosity=k/omega 210! 211 b(i)=b(i)+rhovol*vel(i,6)*(( 212 & (2.d0*(gradvel(1,1,i)**2+gradvel(2,2,i)**2+ 213 & gradvel(3,3,i)**2)+ 214 & (gradvel(1,2,i)+gradvel(2,1,i))**2+ 215 & (gradvel(1,3,i)+gradvel(3,1,i))**2+ 216 & (gradvel(2,3,i)+gradvel(3,2,i))**2))/vel(i,7)) 217 else 218! 219! SST model: other definition of the turbulent 220! viscosity 221! 222 b(i)=b(i)+rhovol*0.31d0*vel(i,6)*(( 223 & (2.d0*(gradvel(1,1,i)**2+gradvel(2,2,i)**2+ 224 & gradvel(3,3,i)**2)+ 225 & (gradvel(1,2,i)+gradvel(2,1,i))**2+ 226 & (gradvel(1,3,i)+gradvel(3,1,i))**2+ 227 & (gradvel(2,3,i)+gradvel(3,2,i))**2))/ 228 & max(0.31d0*vel(i,7),of2(i))) 229 endif 230! 231! transient term 232! 233 constant=rhovol*sc(i) 234 b(i)=b(i)-(a2*velo(i,6)+a3*veloo(i,6))*constant 235 constant=a1*constant 236 ad(i)=ad(i)+constant 237! 238 enddo 239! 240 return 241 end 242