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 e_c3d_v1rhs(co,nk,konl,lakonl,p1,p2,omx,bodyfx, 20 & nbody,ff,nelem,nmethod,rhcon,nrhcon,ielmat,ntmat_,vold,vcon, 21 & dtimef,mi,ttime,time,istep,shcon,nshcon, 22 & iturbulent,nelemface,sideface,nface,compressible, 23 & ipvar,var,ipvarf,varf,ithermal,ipface,nelemload, 24 & sideload,xload,nload,ifreesurface,depth,dgravity,cocon, 25 & ncocon,ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc,iinc, 26 & theta1,bb,physcon,reltimef,xloadold) 27! 28! computation of the velocity element matrix and rhs for the element with 29! element with the topology in konl: step 1 (correction *) 30! 31! ff: rhs 32! 33 implicit none 34! 35 character*1 sideface(*) 36 character*8 lakonl 37 character*20 sideload(*) 38! 39 integer konl(8),ifaceq(8,6),nk,nbody,nelem,ithermal(*),mi(*), 40 & i,j,k,i1,i2,j1,nmethod,ii,jj,id,ipointer,idf, 41 & ig,kk,nrhcon(*),ielmat(mi(3),*),nshcon(*),ntmat_,nope, 42 & nopes,imat,iturbulent,compressible,ipface(*),nelemload(2,*), 43 & mint2d,mint3d,ifacet(6,4),ifacew(8,5),istep,nload, 44 & k1,nelemface(*),nface,ipvar(*),index,ipvarf(*),igl, 45 & ifreesurface,ncocon(2,*),ipompc(*),nodempc(3,*),nmpc, 46 & ikmpc(*),ilmpc(*),iinc,iscale,jltyp,nfield,iemchange, 47 & iflux,node 48! 49 real*8 co(3,*),shp(4,8),dvi,p1(3),p2(3),x2d3,dep,dvel, 50 & bodyfx(3),ff(0:mi(2),8),bf(3),q(3),c1,c2,xsjmod,dtimef2,fric, 51 & rhcon(0:1,ntmat_,*),vel(3),div,shcon(0:3,ntmat_,*),cp, 52 & voldl(0:mi(2),8),xsj2(3),shp2(7,8),omcor,shps(8),ctuf, 53 & vold(0:mi(2),*),om,omx,const,xsj,temp,tt(3),areaj,enthalpy, 54 & vcon(nk,0:mi(2)),vconl(0:mi(2),8),rho,shpv(8),t(3,3), 55 & cvel(3),vkl(3,3),corio(3),xkin,umttot,xload(2,*),f1,f1m, 56 & xtuf,vort,y,f2,unt,a1,arg2,arg1,gamm, 57 & var(*),varf(*),tu,depth(*),dgravity,gamm1, 58 & beta,beta1,beta2,betas,bfv,c3,c4,cdktuf,ckin,cond,gamm2, 59 & cocon(0:6,ntmat_,*),coefmpc(*),press,skin,skin1,skin2,stuf, 60 & stuf1,stuf2,theta1,tuk,turbprandl,tut,umsk,umst,umt,xkappa, 61 & xtu,tvar(2),rhovel(3),dtem(3),dpress(3),aux(3),coords(3), 62 & bb(3,8),tv(3),pgauss(3),physcon(*),dxkin(3),dxtuf(3), 63 & dxsj2,field,reltimef,sinktemp,tvn,xloadold(2,*),tvnk,tvnt, 64 & xsjmodk,xsjmodt,tvk(3),tvt(3) 65! 66 real*8 dtimef,ttime,time 67! 68 ifaceq=reshape((/4,3,2,1,11,10,9,12, 69 & 5,6,7,8,13,14,15,16, 70 & 1,2,6,5,9,18,13,17, 71 & 2,3,7,6,10,19,14,18, 72 & 3,4,8,7,11,20,15,19, 73 & 4,1,5,8,12,17,16,20/),(/8,6/)) 74 ifacet=reshape((/1,3,2,7,6,5, 75 & 1,2,4,5,9,8, 76 & 2,3,4,6,10,9, 77 & 1,4,3,8,10,7/),(/6,4/)) 78 ifacew=reshape((/1,3,2,9,8,7,0,0, 79 & 4,5,6,10,11,12,0,0, 80 & 1,2,5,4,7,14,10,13, 81 & 2,3,6,5,8,15,11,14, 82 & 4,6,3,1,12,15,9,13/),(/8,5/)) 83! 84 tvar(1)=time 85 tvar(2)=ttime+dtimef 86! 87! parameter to express 2.d0/3.d0 88! 89 x2d3=2.d0/3.d0 90 dtimef2=dtimef/2.d0 91! 92! turbulence constants (SST: iturbulent=4) 93! 94 if(iturbulent.gt.0) then 95 a1=0.31d0 96 if(iturbulent.eq.4) then 97 skin1=0.85d0 98 else 99 skin1=0.5d0 100 endif 101 skin2=1.d0 102 stuf1=0.5d0 103 stuf2=0.856d0 104 beta1=0.075d0 105 beta2=0.0828d0 106 betas=0.09d0 107 xkappa= 0.41d0 108! 109 gamm1=beta1/betas-stuf1*xkappa*xkappa/dsqrt(betas) 110 gamm2=beta2/betas-stuf2*xkappa*xkappa/dsqrt(betas) 111! 112 xtu=10.d0*physcon(5)/physcon(8) 113 xtu=xtu*xtu 114 c3=betas*xtu*10.d0**(-3.5d0) 115 endif 116! 117 imat=ielmat(1,nelem) 118! 119 if(lakonl(4:4).eq.'4') then 120 nope=4 121 mint3d=1 122 elseif(lakonl(4:4).eq.'6') then 123 nope=6 124 mint3d=2 125 elseif(lakonl(4:5).eq.'8R') then 126 nope=8 127 mint3d=1 128 elseif(lakonl(4:4).eq.'8') then 129 nope=8 130 mint3d=8 131 endif 132! 133! initialisation for distributed forces 134! 135 do i1=1,nope 136 do i2=0,mi(2) 137 ff(i2,i1)=0.d0 138 enddo 139 enddo 140 do i1=1,nope 141 do i2=1,3 142 bb(i2,i1)=0.d0 143 enddo 144 enddo 145! 146! temperature, velocity, conservative variables 147! (rho*velocity and rho) and if turbulent 148! rho*turbulence variables 149! 150 do i1=1,nope 151 do i2=0,mi(2) 152 voldl(i2,i1)=vold(i2,konl(i1)) 153 enddo 154 do i2=0,mi(2) 155 vconl(i2,i1)=vcon(konl(i1),i2) 156 enddo 157c if(nelem.eq.1) then 158c write(*,*) 'voldl(5..6) ',voldl(5,i1),voldl(6,i1) 159c write(*,*) 'vconl(5..6) ',vconl(5,i1),vconl(6,i1) 160c endif 161 enddo 162! 163! computation of the matrix: loop over the Gauss points 164! 165 index=ipvar(nelem) 166 do kk=1,mint3d 167! 168! copying the shape functions, their derivatives and the 169! Jacobian determinant from field var 170! 171 do jj=1,nope 172 do ii=1,4 173 index=index+1 174 shp(ii,jj)=var(index) 175 enddo 176 enddo 177 index=index+1 178 xsj=var(index) 179 index=index+1 180 y=var(index) 181! 182 xsjmod=dtimef*xsj 183! 184! the temperature temp 185! the velocity vel(*) 186! rho times the velocity rhovel(*) 187! the temperature gradient dtem(*) 188! the velocity gradient vkl(*,*) 189! the pressure gradient dpress(*) 190! 191 temp=0.d0 192 do j1=1,3 193 vel(j1)=0.d0 194 rhovel(j1)=0.d0 195 dtem(j1)=0.d0 196 do k1=1,3 197 vkl(j1,k1)=0.d0 198 enddo 199 dpress(j1)=0.d0 200 enddo 201! 202 do i1=1,nope 203 temp=temp+shp(4,i1)*voldl(0,i1) 204 do j1=1,3 205 vel(j1)=vel(j1)+shp(4,i1)*voldl(j1,i1) 206 rhovel(j1)=rhovel(j1)+shp(4,i1)*vconl(j1,i1) 207 dtem(j1)=dtem(j1)+shp(j1,i1)*voldl(0,i1) 208 do k1=1,3 209 vkl(j1,k1)=vkl(j1,k1)+shp(k1,i1)*voldl(j1,i1) 210 enddo 211 dpress(j1)=dpress(j1)+shp(j1,i1)*voldl(4,i1) 212 enddo 213 enddo 214! 215! the divergence of the velocity div 216! 217 if(compressible.eq.1) then 218 div=vkl(1,1)+vkl(2,2)+vkl(3,3) 219 else 220 div=0.d0 221 endif 222! 223! the divergence of the shape function times the velocity shpv(*) 224! the convective enthalpy 225! the convective velocity cvel(*) 226! 227 shpv(1)=0.d0 228 enthalpy=0.d0 229 do j1=1,3 230 cvel(j1)=0.d0 231 enddo 232! 233 do i1=1,nope 234 shpv(i1)=shp(1,i1)*vel(1)+shp(2,i1)*vel(2)+ 235 & shp(3,i1)*vel(3)+shp(4,i1)*div 236 enthalpy=enthalpy+shpv(i1)*(vconl(0,i1)+voldl(4,i1)) 237 do j1=1,3 238 cvel(j1)=cvel(j1)+shpv(i1)*vconl(j1,i1) 239 enddo 240 enddo 241! 242! creating auxiliary variable shps 243! 244 do i1=1,nope 245 shps(i1)=xsjmod*(shp(4,i1)+dtimef2*shpv(i1)) 246 enddo 247! 248! material data (dynamic viscosity) 249! 250 call materialdata_dvifem(imat,ntmat_,temp,shcon,nshcon,dvi) 251! 252! determining the dissipative stress 253! 254 do i1=1,3 255 do j1=i1,3 256 t(i1,j1)=vkl(i1,j1)+vkl(j1,i1) 257 enddo 258 if(compressible.eq.1) t(i1,i1)=t(i1,i1)-x2d3*div 259 enddo 260! 261! calculation of the density in case of body forces and/or 262! turbulence 263! 264 if((nbody.ne.0).or.(iturbulent.ne.0)) then 265 if(compressible.eq.1) then 266! 267! gas 268! 269 rho=0.d0 270 do i1=1,nope 271 rho=rho+shp(4,i1)*vconl(4,i1) 272 enddo 273 else 274! 275! liquid 276! 277 call materialdata_rho(rhcon,nrhcon,imat,rho, 278 & temp,ntmat_,ithermal) 279 endif 280 endif 281! 282! calculation of the turbulent kinetic energy, turbulence 283! frequency and their spatial derivatives for gases and liquids 284! 285 if(iturbulent.gt.0) then 286 xkin=0.d0 287 xtuf=0.d0 288 do i1=1,nope 289 xkin=xkin+shp(4,i1)*voldl(5,i1) 290 xtuf=xtuf+shp(4,i1)*voldl(6,i1) 291 enddo 292! 293! adding the turbulent stress 294! 295! factor F2 296! 297 c1=dsqrt(xkin)/(0.09d0*xtuf*y) 298 c2=500.d0*dvi/(y*y*xtuf*rho) 299! 300! kinematic turbulent viscosity 301! 302 if(iturbulent.eq.4) then 303! 304! vorticity 305! 306 vort=dsqrt((vkl(3,2)-vkl(2,3))**2+ 307 & (vkl(1,3)-vkl(3,1))**2+ 308 & (vkl(2,1)-vkl(1,2))**2) 309 arg2=max(2.d0*c1,c2) 310 f2=dtanh(arg2*arg2) 311 unt=a1*xkin/max(a1*xtuf,vort*f2) 312 else 313 unt=xkin/xtuf 314 endif 315! 316! calculating the production (anisotropic part of 317! the turbulent stress is, apart from the dynamic 318! viscosity, identical to the viscous stress) 319! 320 tu=t(1,1)*vkl(1,1)+t(1,2)*t(1,2)+t(1,3)*t(1,3)+ 321 & t(2,2)*vkl(2,2)+t(2,3)*t(2,3)+t(3,3)*vkl(3,3) 322! 323! correction for compressible fluids 324! 325 if(compressible.eq.1) then 326 tu=tu-x2d3*xkin*div/unt 327 endif 328! 329! calculating the turbulent stress 330! 331 umttot=dvi+unt*rho 332 do i1=1,3 333 do j1=i1,3 334 t(i1,j1)=umttot*t(i1,j1) 335 enddo 336 t(i1,i1)=t(i1,i1)-x2d3*rho*xkin 337 enddo 338 else 339! 340 do i1=1,3 341 do j1=i1,3 342 t(i1,j1)=dvi*t(i1,j1) 343 enddo 344 enddo 345 endif 346! 347! 1a. rhs of the first part of the momentum equation 348! 349 do jj=1,nope 350! 351! convective + diffusive 352! 353 ff(1,jj)=ff(1,jj)-cvel(1)*shps(jj)-xsjmod* 354 & (shp(1,jj)*t(1,1)+shp(2,jj)*t(1,2)+shp(3,jj)*t(1,3)) 355 ff(2,jj)=ff(2,jj)-cvel(2)*shps(jj)-xsjmod* 356 & (shp(1,jj)*t(1,2)+shp(2,jj)*t(2,2)+shp(3,jj)*t(2,3)) 357 ff(3,jj)=ff(3,jj)-cvel(3)*shps(jj)-xsjmod* 358 & (shp(1,jj)*t(1,3)+shp(2,jj)*t(2,3)+shp(3,jj)*t(3,3)) 359 enddo 360! 361! computation of contribution due to body forces 362! 363 if(nbody.ne.0) then 364! 365! initialisation for the body forces 366! 367 om=omx*rho 368 omcor=2.d0*rho*dsqrt(omx) 369! 370 if(om.gt.0.d0) then 371 do i1=1,3 372! 373! computation of the global coordinates of the gauss 374! point 375! 376 q(i1)=0.d0 377 do j1=1,nope 378 q(i1)=q(i1)+shp(4,j1)*co(i1,konl(j1)) 379 enddo 380! 381 q(i1)=q(i1)-p1(i1) 382 enddo 383 const=q(1)*p2(1)+q(2)*p2(2)+q(3)*p2(3) 384! 385! Coriolis forces 386! 387 omcor=2.d0*rho*dsqrt(omx) 388 corio(1)=vel(2)*p2(3)-vel(3)*p2(2) 389 corio(2)=vel(3)*p2(1)-vel(1)*p2(3) 390 corio(3)=vel(1)*p2(2)-vel(2)*p2(1) 391 endif 392! 393 if(ifreesurface.eq.0) then 394 do ii=1,3 395 bf(ii)=bodyfx(ii)*rho 396 enddo 397! 398! inclusion of the centrifugal force into the body force 399! 400 if(om.gt.0.d0) then 401 do i1=1,3 402 bf(i1)=bf(i1)+(q(i1)-const*p2(i1))*om+ 403 & corio(i1)*omcor 404 enddo 405 endif 406 else 407! 408! shallow water calculation 409! effect of varying depth; 410! the effect of the centrifugal force on dgravity is neglected 411! 412 dep=0.d0 413 do j1=1,3 414 bf(j1)=0.d0 415 enddo 416! 417 do i1=1,nope 418 dep=dep+shp(4,i1)*depth(konl(i1)) 419 do j1=1,3 420 bf(j1)=bf(j1)+shp(j1,i1)*depth(konl(i1)) 421 enddo 422 enddo 423 do j1=1,3 424 bf(j1)=bf(j1)*(rho-dep)*dgravity 425 enddo 426! 427 if(om.gt.0.d0) then 428 do i1=1,2 429 bf(i1)=bf(i1)+(q(i1)-const*p2(i1))*om+ 430 & corio(i1)*omcor 431 enddo 432 endif 433! 434! bottom friction 435! 436c fric=0.02d0 437 fric=0.01d0 438 dvel=dsqrt(vel(1)*vel(1)+vel(2)*vel(2)+vel(3)*vel(3)) 439 do j1=1,3 440 bf(j1)=bf(j1)-fric*dvel*vel(j1)/8.d0 441 enddo 442 endif 443! 444! storing the body force 445! 446 bfv=bf(1)*vel(1)+bf(2)*vel(2)+bf(3)*vel(3) 447! 448! 1b. rhs of the first part of the momentum equation: 449! body force contribution 450! 451 do jj=1,nope 452 ff(1,jj)=ff(1,jj)+bf(1)*shps(jj) 453 ff(2,jj)=ff(2,jj)+bf(2)*shps(jj) 454 ff(3,jj)=ff(3,jj)+bf(3)*shps(jj) 455 enddo 456 endif 457! 458! 2. rhs of the mass equation 459! 460 do j1=1,3 461 aux(j1)=xsjmod*(rhovel(j1)-dtimef*theta1*dpress(j1)) 462 enddo 463! 464 do jj=1,nope 465 ff(4,jj)=ff(4,jj)+ 466 & shp(1,jj)*aux(1)+shp(2,jj)*aux(2)+shp(3,jj)*aux(3) 467 enddo 468! 469! 3. rhs of the second part of the momentum equation: 470! 471 if(compressible.eq.1) then 472! 473! explicit compressible 474! 475 do jj=1,nope 476 bb(1,jj)=bb(1,jj)-dpress(1)*shps(jj) 477 bb(2,jj)=bb(2,jj)-dpress(2)*shps(jj) 478 bb(3,jj)=bb(3,jj)-dpress(3)*shps(jj) 479 enddo 480 else 481! 482! implicit incompressible 483! 484 do jj=1,nope 485 bb(1,jj)=bb(1,jj)-xsjmod*shp(4,jj)*dpress(1) 486 bb(2,jj)=bb(2,jj)-xsjmod*shp(4,jj)*dpress(2) 487 bb(3,jj)=bb(3,jj)-xsjmod*shp(4,jj)*dpress(3) 488 enddo 489 endif 490! 491! 4. rhs of the energy equation: 492! 493 if(ithermal(1).gt.0) then 494! 495! viscous conductivity 496! 497 call materialdata_cond(imat,ntmat_,temp,cocon,ncocon,cond) 498! 499! adding the turbulent conductivity 500! 501 if(iturbulent.gt.0) then 502 call materialdata_cp(imat,ntmat_,temp,shcon,nshcon,cp) 503 turbprandl=0.9d0 504 cond=cond+cp*rho*unt/turbprandl 505 endif 506! 507! calculating the total dissipative stress x velocity 508! (viscous + turbulent) 509! 510 tv(1)=t(1,1)*vel(1)+t(1,2)*vel(2)+t(1,3)*vel(3) 511 tv(2)=t(1,2)*vel(1)+t(2,2)*vel(2)+t(2,3)*vel(3) 512 tv(3)=t(1,3)*vel(1)+t(2,3)*vel(2)+t(3,3)*vel(3) 513! 514! determining stress x velocity + conductivity x 515! temperature gradient 516! 517 do i1=1,3 518 tv(i1)=tv(i1)+cond*dtem(i1) 519 enddo 520! 521! determination of the rhs of the energy equations 522! 523 do jj=1,nope 524 ff(0,jj)=ff(0,jj)-shps(jj)*enthalpy-xsjmod* 525 & (shp(1,jj)*tv(1)+shp(2,jj)*tv(2)+shp(3,jj)*tv(3)) 526 enddo 527! 528! computation of contribution due to body forces 529! 530 if(nbody.ne.0) then 531 do jj=1,nope 532 ff(0,jj)=ff(0,jj)+shps(jj)*bfv 533 enddo 534 endif 535! 536! distributed heat flux 537! 538 if(nload.gt.0) then 539 call nident2(nelemload,nelem,nload,id) 540 areaj=xsj 541 do 542 if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit 543 if(sideload(id)(1:2).ne.'BF') then 544 id=id-1 545 cycle 546 endif 547 if(sideload(id)(3:4).eq.'NU') then 548 do j=1,3 549 pgauss(j)=0.d0 550 do i1=1,nope 551 pgauss(j)=pgauss(j)+ 552 & shp(4,i1)*co(j,konl(i1)) 553 enddo 554 enddo 555 jltyp=1 556 iscale=1 557 call dflux(xload(1,id),temp,istep,iinc,tvar, 558 & nelem,kk,pgauss,jltyp,temp,press,sideload(id), 559 & areaj,vold,co,lakonl,konl,ipompc,nodempc,coefmpc, 560 & nmpc,ikmpc,ilmpc,iscale,mi) 561 endif 562 do jj=1,nope 563 ff(0,jj)=ff(0,jj)+shps(jj)*xload(1,id) 564 enddo 565 exit 566 enddo 567 endif 568 endif 569! 570! 5. rhs of the turbulence equations: 571! 572 if(iturbulent.gt.0) then 573! 574! convective turbulent kinetic energy: ckin 575! convective turbulence frequency: ctuf 576! 577 ckin=0.d0 578 ctuf=0.d0 579 do i1=1,nope 580 ckin=ckin+shpv(i1)*vconl(5,i1) 581 ctuf=ctuf+shpv(i1)*vconl(6,i1) 582 enddo 583c if(nelem.eq.1) then 584c write(*,*) 'ckin ',ckin 585c write(*,*) 'ctuf ',ctuf 586c endif 587! 588! gradient of k and omega 589! 590 do j1=1,3 591 dxkin(j1)=0.d0 592 dxtuf(j1)=0.d0 593 enddo 594 do i1=1,nope 595 do j1=1,3 596 dxkin(j1)=dxkin(j1)+shp(j1,i1)*voldl(5,i1) 597 dxtuf(j1)=dxtuf(j1)+shp(j1,i1)*voldl(6,i1) 598 enddo 599 enddo 600c if(nelem.eq.1) then 601c write(*,*) 'dxkin ',(dxkin(j1),j1=1,3) 602c write(*,*) 'dxtuf ',(dxtuf(j1),j1=1,3) 603c endif 604! 605! auxiliary variable 606! 607 c4=2.d0*rho*stuf2* 608 & (dxkin(1)*dxtuf(1)+dxkin(2)*dxtuf(2)+dxkin(3)*dxtuf(3))/ 609 & xtuf 610! 611! dynamic turbulent viscosity 612! 613 umt=unt*rho 614! 615! factor F1 616! 617 if(iturbulent.eq.1) then 618! 619! k-epsilon model 620! 621 f1=0.d0 622 elseif(iturbulent.eq.2) then 623! 624! k-omega model 625! 626 f1=1.d0 627 else 628! 629! BSL/SST model 630! 631 cdktuf=max(c4,1.d-20) 632 arg1=min(max(c1,c2),4.d0*rho*stuf2*xkin/(cdktuf*y*y)) 633 f1=dtanh(arg1**4.d0) 634 endif 635 f1m=1.d0-f1 636! 637! interpolation of the constants 638! 639 skin=f1*skin1+f1m*skin2 640 stuf=f1*stuf1+f1m*stuf2 641 beta=f1*beta1+f1m*beta2 642 gamm=f1*gamm1+f1m*gamm2 643! 644! source terms: productivity - dissipation 645! 646 umsk=dvi+skin*umt 647 umst=dvi+stuf*umt 648! 649! production limiter active: P=unt*tu<=20*betas*k*omega 650! Menter, F.R., "Zonal Two Equation k-omega Turbulence Models for 651! Aerodynamic Flows," AIAA Paper 93-2906, July 1993. 652! 653 tuk=rho*(unt*tu-betas*xtuf*xkin) 654 tut=rho*(gamm*tu-beta*xtuf*xtuf)+f1m*c4 655! 656! add controlled decay 657! Spalart, P.R. and Rumsey, C.L., "Effective Inflow Conditions for 658! Turbulence Models in Aerodynamic Calculations," AIAA Journal, 659! Vol. 45, No. 10, 2007,pp.2544-2553. 660! 661 tuk=tuk+c3*dvi 662 tut=tut+beta*xtu*rho 663! 664 do i1=1,3 665 dxkin(i1)=dxkin(i1)*umsk 666 dxtuf(i1)=dxtuf(i1)*umst 667 enddo 668! 669! determination of rhs 670! 671 do jj=1,nope 672! 673 ff(5,jj)=ff(5,jj)-shps(jj)*(ckin-tuk)-xsjmod* 674 & (shp(1,jj)*dxkin(1)+shp(2,jj)*dxkin(2) 675 & +shp(3,jj)*dxkin(3)) 676 ff(6,jj)=ff(6,jj)-shps(jj)*(ctuf-tut)-xsjmod* 677 & (shp(1,jj)*dxtuf(1)+shp(2,jj)*dxtuf(2) 678 & +shp(3,jj)*dxtuf(3)) 679 enddo 680 endif 681! 682 enddo 683! 684! area integrals 685! 686 if(nface.ne.0) then 687 index=ipvarf(nelem) 688c write(*,*) 'e_c3d_v1rhs nelem ipvarf(nelem)',nelem,index 689! 690! external boundaries 691! 692 nopes=0 693 idf=ipface(nelem) 694 do 695 if((idf.eq.0).or.(nelemface(idf).ne.nelem)) exit 696 ig=ichar(sideface(idf)(1:1))-48 697! 698! check for distributed flux 699! an adiabatic face must be declared as a face with 700! distributed flux zero! 701! 702 iflux=0 703 call nident2(nelemload,nelem,nload,id) 704 do 705 if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit 706 if((sideload(id)(1:1).ne.'F').and. 707 & (sideload(id)(1:1).ne.'R').and. 708 & (sideload(id)(1:1).ne.'S')) then 709 id=id-1 710 cycle 711 endif 712 igl=ichar(sideload(id)(2:2))-48 713 if(igl.ne.ig) then 714 id=id-1 715 cycle 716 endif 717 iflux=1 718 exit 719 enddo 720! 721 if(nopes.eq.0) then 722 if(lakonl(4:4).eq.'4') then 723 nopes=3 724 mint2d=1 725 elseif(lakonl(4:4).eq.'6') then 726 mint2d=1 727 elseif(lakonl(4:5).eq.'8R') then 728 nopes=4 729 mint2d=1 730 elseif(lakonl(4:4).eq.'8') then 731 nopes=4 732 mint2d=4 733 endif 734 endif 735! 736 if(lakonl(4:4).eq.'6') then 737 if(ig.le.2) then 738 nopes=3 739 else 740 nopes=4 741 endif 742 endif 743! 744c write(*,*) 'e_c3d_v1rhs ',index,4*nope+nopes+4 745 do i=1,mint2d 746! 747! facial shape functions 748! local surface normal 749! 750 do i1=1,nopes 751 index=index+1 752 shp2(4,i1)=varf(index) 753 enddo 754 do i1=1,3 755 index=index+1 756 xsj2(i1)=varf(index) 757 enddo 758! 759! derivative of the volumetric shape functions 760! needed for the temperature, velocity gradients and 761! gradients of k and omega (turbulence) 762! 763 do i1=1,nope 764 do j1=1,4 765 index=index+1 766 shp(j1,i1)=varf(index) 767 enddo 768 enddo 769 index=index+1 770 y=varf(index) 771! 772! calculating of 773! the temperature temp 774! the velocity vel(*) 775! rho times the velocity rhovel(*) 776! the velocity gradient vkl 777! 778 temp=0.d0 779 do j1=1,3 780 vel(j1)=0.d0 781 rhovel(j1)=0.d0 782 do k1=1,3 783 vkl(j1,k1)=0.d0 784 enddo 785 enddo 786! 787 do i1=1,nope 788 temp=temp+shp(4,i1)*voldl(0,i1) 789 do j1=1,3 790 vel(j1)=vel(j1)+shp(4,i1)*voldl(j1,i1) 791c write(*,*) 'e_c3d_v1rhs ',shp(4,i1) 792 rhovel(j1)=rhovel(j1)+shp(4,i1)*vconl(j1,i1) 793 do k1=1,3 794 vkl(j1,k1)=vkl(j1,k1)+shp(k1,i1)*voldl(j1,i1) 795 enddo 796 enddo 797 enddo 798! 799 if(iflux.eq.0) then 800! 801! calculating of the temperature gradient dtem 802! in the integration point 803! 804 do j1=1,3 805 dtem(j1)=0.d0 806 enddo 807 do i1=1,nope 808 do j1=1,3 809 dtem(j1)=dtem(j1)+shp(j1,i1)*voldl(0,i1) 810 enddo 811 enddo 812 endif 813! 814 if(compressible.eq.1) then 815 div=vkl(1,1)+vkl(2,2)+vkl(3,3) 816 else 817 div=0.d0 818 endif 819! 820! material data (dynamic viscosity) 821! 822 call materialdata_dvifem(imat,ntmat_,temp,shcon,nshcon, 823 & dvi) 824! 825! determining the dissipative stress 826! 827 do i1=1,3 828 do j1=i1,3 829 t(i1,j1)=vkl(i1,j1)+vkl(j1,i1) 830 enddo 831 if(compressible.eq.1) t(i1,i1)=t(i1,i1)-x2d3*div 832 enddo 833! 834! calculation of the density for gases 835! 836! calculation of the turbulent kinetic energy, turbulence 837! frequency and their spatial derivatives for gases and liquids 838! 839 if(iturbulent.gt.0) then 840 if(compressible.eq.1) then 841! 842! gas 843! 844 rho=0.d0 845 do i1=1,nope 846 rho=rho+shp(4,i1)*vconl(4,i1) 847 enddo 848 else 849! 850! liquid 851! 852 call materialdata_rho(rhcon,nrhcon,imat,rho, 853 & temp,ntmat_,ithermal) 854! 855! calculation of k, omega an y 856! 857 endif 858 xkin=0.d0 859 xtuf=0.d0 860 do i1=1,nope 861 xkin=xkin+shp(4,i1)*voldl(5,i1) 862 xtuf=xtuf+shp(4,i1)*voldl(6,i1) 863 enddo 864! 865! calculation of turbulent auxiliary variables 866! 867! factor F2 868! 869 if(y.gt.0.d0) then 870 c1=dsqrt(xkin)/(0.09d0*xtuf*y) 871 c2=500.d0*dvi/(y*y*xtuf*rho) 872 endif 873! 874! kinematic and dynamic turbulent viscosity 875! 876 if(iturbulent.eq.4) then 877! 878! vorticity 879! 880 vort=dsqrt((vkl(3,2)-vkl(2,3))**2+ 881 & (vkl(1,3)-vkl(3,1))**2+ 882 & (vkl(2,1)-vkl(1,2))**2) 883 if(y.gt.0.d0) then 884 arg2=max(2.d0*c1,c2) 885 f2=dtanh(arg2*arg2) 886 else 887 f2=1.d0 888 endif 889 unt=a1*xkin/max(a1*xtuf,vort*f2) 890 else 891 unt=xkin/xtuf 892 endif 893! 894 umttot=dvi+unt*rho 895 do i1=1,3 896 do j1=i1,3 897 t(i1,j1)=umttot*t(i1,j1) 898 enddo 899 t(i1,i1)=t(i1,i1)-x2d3*rho*xkin 900 enddo 901 else 902! 903 do i1=1,3 904 do j1=i1,3 905 t(i1,j1)=dvi*t(i1,j1) 906 enddo 907 enddo 908 endif 909! 910! stress vector 911! 912 tt(1)=(t(1,1)*xsj2(1)+t(1,2)*xsj2(2)+t(1,3)*xsj2(3))* 913 & dtimef 914 tt(2)=(t(1,2)*xsj2(1)+t(2,2)*xsj2(2)+t(2,3)*xsj2(3))* 915 & dtimef 916 tt(3)=(t(1,3)*xsj2(1)+t(2,3)*xsj2(2)+t(3,3)*xsj2(3))* 917 & dtimef 918! 919! stress x velocity 920! 921 tv(1)=t(1,1)*vel(1)+t(1,2)*vel(2)+t(1,3)*vel(3) 922 tv(2)=t(1,2)*vel(1)+t(2,2)*vel(2)+t(2,3)*vel(3) 923 tv(3)=t(1,3)*vel(1)+t(2,3)*vel(2)+t(3,3)*vel(3) 924! 925! adding conductivity in case the flux is not given by a 926! *DFLUX, *FILM or *RADIATE card 927! 928 if(iflux.eq.0) then 929 do i1=1,3 930 tv(i1)=tv(i1)+cond*dtem(i1) 931 enddo 932 endif 933! 934 tvn=tv(1)*xsj2(1)+tv(2)*xsj2(2)+tv(3)*xsj2(3) 935! 936! modifying tvn in case of a *DFLUX, *FILM or *RADIATE 937! card 938! 939 if(iflux.eq.1) then 940 dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 941 & xsj2(3)*xsj2(3)) 942 areaj=dxsj2 943 sinktemp=xload(2,id) 944! 945! for nonuniform load: determine the coordinates of the 946! point (transferred into the user subroutine) 947! 948 if((sideload(id)(3:4).eq.'NU').or. 949 & (sideload(id)(5:6).eq.'NU')) then 950 if(nope.eq.8) then 951 do k=1,3 952 coords(k)=0.d0 953 do j=1,nopes 954 coords(k)=coords(k)+ 955 & co(k,konl(ifaceq(j,ig)))*shp2(4,j) 956 enddo 957 enddo 958 elseif(nope.eq.4) then 959 do k=1,3 960 coords(k)=0.d0 961 do j=1,nopes 962 coords(k)=coords(k)+ 963 & co(k,konl(ifacet(j,ig)))*shp2(4,j) 964 enddo 965 enddo 966 else 967 do k=1,3 968 coords(k)=0.d0 969 do j=1,nopes 970 coords(k)=coords(k)+ 971 & co(k,konl(ifacew(j,ig)))*shp2(4,j) 972 enddo 973 enddo 974 endif 975 jltyp=ichar(sideload(id)(2:2))-48 976 jltyp=jltyp+10 977 if(sideload(id)(1:1).eq.'S') then 978 iscale=1 979 call dflux(xload(1,id),temp,istep,iinc,tvar, 980 & nelem,i,coords,jltyp,temp,press, 981 & sideload(id),areaj,vold,co,lakonl,konl, 982 & ipompc,nodempc,coefmpc,nmpc,ikmpc,ilmpc, 983 & iscale,mi) 984 if((nmethod.eq.1).and.(iscale.ne.0)) 985 & xload(1,id)=xloadold(1,id)+ 986 & (xload(1,id)-xloadold(1,id))*reltimef 987 elseif(sideload(id)(1:1).eq.'F') then 988 call film(xload(1,id),sinktemp,temp,istep, 989 & iinc,tvar,nelem,i,coords,jltyp,field, 990 & nfield,sideload(id),node,areaj,vold,mi) 991 if(nmethod.eq.1) xload(1,id)=xloadold(1,id)+ 992 & (xload(1,id)-xloadold(1,id))*reltimef 993 elseif(sideload(id)(1:1).eq.'R') then 994 call radiate(xload(1,id),xload(2,id),temp,istep, 995 & iinc,tvar,nelem,i,coords,jltyp,field, 996 & nfield,sideload(id),node,areaj,vold,mi, 997 & iemchange) 998 if(nmethod.eq.1) xload(1,id)=xloadold(1,id)+ 999 & (xload(1,id)-xloadold(1,id))*reltimef 1000 endif 1001 endif 1002! 1003 if(sideload(id)(1:1).eq.'S') then 1004! 1005! flux INTO the face is positive (input deck convention) 1006! this is different from the convention in the theory 1007! 1008 tvn=tvn+xload(1,id)*dxsj2 1009 elseif(sideload(id)(1:1).eq.'F') then 1010 tvn=tvn-xload(1,id)*(temp-sinktemp)*dxsj2 1011 elseif(sideload(id)(1:1).eq.'R') then 1012 tvn=tvn-physcon(2)* 1013 & xload(1,id)*((temp-physcon(1))**4- 1014 & (xload(2,id)-physcon(1))**4)*dxsj2 1015 endif 1016 endif 1017! 1018 xsjmod=tvn*dtimef 1019! 1020 if(iturbulent.gt.0) then 1021! 1022! calculation of the spatial derivatives of the turbulent kinetic energy 1023! and the turbulence frequency for gases and liquids 1024! 1025 do j1=1,3 1026 dxkin(j1)=0.d0 1027 dxtuf(j1)=0.d0 1028 enddo 1029 do i1=1,nope 1030 do j1=1,3 1031 dxkin(j1)=dxkin(j1)+shp(j1,i1)*voldl(5,i1) 1032 dxtuf(j1)=dxtuf(j1)+shp(j1,i1)*voldl(6,i1) 1033 enddo 1034 enddo 1035! 1036! auxiliary variable 1037! 1038 c4=2.d0*rho*stuf2* 1039 & (dxkin(1)*dxtuf(1)+dxkin(2)*dxtuf(2)+ 1040 & dxkin(3)*dxtuf(3))/xtuf 1041! 1042! dynamic turbulent viscosity 1043! 1044 umt=unt*rho 1045! 1046! factor F1 1047! 1048 if(iturbulent.eq.1) then 1049! 1050! k-epsilon model 1051! 1052 f1=0.d0 1053 elseif(iturbulent.eq.2) then 1054! 1055! k-omega model 1056! 1057 f1=1.d0 1058 else 1059! 1060! BSL/SST model 1061! 1062 if(y.gt.0.d0) then 1063! 1064! finite distance from wall 1065! 1066 cdktuf=max(c4,1.d-20) 1067 arg1= 1068 & min(max(c1,c2),4.d0*rho*stuf2*xkin/(cdktuf*y*y)) 1069 f1=dtanh(arg1**4.d0) 1070 else 1071! 1072! wall 1073! 1074 f1=1.d0 1075 endif 1076 endif 1077 f1m=1.d0-f1 1078! 1079! interpolation of the constants 1080! 1081 skin=f1*skin1+f1m*skin2 1082 stuf=f1*stuf1+f1m*stuf2 1083! 1084! auxiliary quantities 1085! 1086 umsk=dvi+skin*umt 1087 umst=dvi+stuf*umt 1088! 1089! determining the stress and and stress x velocity + conductivity x 1090! temperature gradient 1091! 1092 do i1=1,3 1093 tvk(i1)=umsk*dxkin(i1) 1094 tvt(i1)=umsk*dxtuf(i1) 1095 enddo 1096! 1097 tvnk=tvk(1)*xsj2(1)+tvk(2)*xsj2(2)+tvk(3)*xsj2(3) 1098 tvnt=tvt(1)*xsj2(1)+tvt(2)*xsj2(2)+tvt(3)*xsj2(3) 1099! 1100 xsjmodk=tvnk*dtimef 1101 xsjmodt=tvnt*dtimef 1102 endif 1103! 1104 do k=1,nopes 1105 if(nope.eq.8) then 1106 ipointer=ifaceq(k,ig) 1107 elseif(nope.eq.4) then 1108 ipointer=ifacet(k,ig) 1109 else 1110 ipointer=ifacew(k,ig) 1111 endif 1112! 1113! 1a. rhs of the first part of the momentum equation 1114! 1115 ff(1,ipointer)=ff(1,ipointer)+shp2(4,k)*tt(1) 1116 ff(2,ipointer)=ff(2,ipointer)+shp2(4,k)*tt(2) 1117 ff(3,ipointer)=ff(3,ipointer)+shp2(4,k)*tt(3) 1118! 1119! 2. rhs of the mass equation 1120! 1121 ff(4,ipointer)=ff(4,ipointer)-shp2(4,k)* 1122 & (rhovel(1)*xsj2(1)+rhovel(2)*xsj2(2)+ 1123 & rhovel(3)*xsj2(3))*dtimef 1124! 1125! 4. rhs of the energy equation: 1126! 1127 if(ithermal(1).gt.0) then 1128 ff(0,ipointer)=ff(0,ipointer)+shp2(4,k)*xsjmod 1129 endif 1130! 1131! 5. rhs of the turbulence equations: 1132! 1133 if(iturbulent.gt.0) then 1134 ff(5,ipointer)=ff(5,ipointer)+shp2(4,k)*xsjmodk 1135 ff(6,ipointer)=ff(6,ipointer)+shp2(4,k)*xsjmodt 1136 endif 1137 enddo 1138 enddo 1139! 1140 idf=idf-1 1141 enddo 1142 endif 1143! 1144c if(nelem.eq.1) then 1145c write(*,*) 'e_c3d_v1rhs' 1146c do k=1,8 1147c write(*,*) nelem,k,(ff(j,k),j=5,6) 1148c enddo 1149c endif 1150 return 1151 end 1152