1 subroutine e3b (ul, yl, acl, iBCB, BCB, 2 & shpb, shglb, 3 & xlb, rl, sgn, dwl, xKebe) 4c 5c---------------------------------------------------------------------- 6c 7c This routine calculates the 3D RHS residual of the fluid boundary 8c elements. 9c 10c input: 11c yl (npro,nshl,ndof) : Y variables 12c iBCB (npro,ndiBCB) : boundary condition code (iBCB(:,1) is 13c a bit tested boundary integral flag i.e. 14c if set to value of BCB if set to floating value 15c iBCB(:,1) : convective flux * 1 0 (ditto to all below) 16c pressure flux * 2 17c viscous flux * 4 18c heat flux * 8 19c turbulence wall * 16 20c scalarI flux * 16*2^I 21c (where I is the scalar number) 22c 23c iBCB(:,2) is the srfID given by the user in MGI that we will 24c collect integrated fluxes for. 25c 26c BCB (npro,nshlb,ndBCB) : Boundary Condition values 27c BCB (1) : mass flux 28c BCB (2) : pressure 29c BCB (3) : viscous flux in x1-direc. 30c BCB (4) : viscous flux in x2-direc. 31c BCB (5) : viscous flux in x3-direc. 32c BCB (6) : heat flux 33c shpb (nen,ngaussb) : boundary element shape-functions 34c shglb (nsd,nen,ngaussb) : boundary element grad-shape-functions 35c xlb (npro,nenl,nsd) : nodal coordinates at current step 36c 37c output: 38c rl (npro,nshl,nflow) : element residual 39c 40c Note: Always the first side of the element is on the boundary. 41c However, note that for higher-order elements the nodes on 42c the boundary side are not the first nshlb nodes, see the 43c array mnodeb. 44c 45c 46c Zdenek Johan, Summer 1990. (Modified from e2b.f) 47c Zdenek Johan, Winter 1991. (Fortran 90) 48c Alberto Figueroa, Winter 2004. CMM-FSI 49c Irene Vignon, Spring 2004 50c---------------------------------------------------------------------- 51c 52 include "common.h" 53c 54 dimension yl(npro,nshl,ndof), iBCB(npro,ndiBCB), 55 & BCB(npro,nshlb,ndBCB), shpb(nshl,ngaussb), 56 & shglb(nsd,nshl,ngaussb), 57 & xlb(npro,nenl,nsd), ul(npro,nshl,nsd), 58 & acl(npro,nshl,ndof), 59 & rl(npro,nshl,nflow) 60c 61 dimension g1yi(npro,ndof), g2yi(npro,ndof), 62 & g3yi(npro,ndof), WdetJb(npro), 63 & bnorm(npro,nsd) 64c 65 dimension u1(npro), u2(npro), 66 & u3(npro), rho(npro), 67 & unm(npro), pres(npro), 68 & vdot(npro,nsd), rlKwall(npro,nshlb,nsd) 69c 70 dimension rou(npro), rmu(npro) 71c 72 dimension tau1n(npro), 73 & tau2n(npro), tau3n(npro) 74c 75 dimension lnode(27), sgn(npro,nshl), 76 & shape(npro,nshl), shdrv(npro,nsd,nshl), 77 & rNa(npro,4) 78 79 real*8 xmudmi(npro,ngauss), dwl(npro,nshl) 80c 81 dimension xKebe(npro,9,nshl,nshl), rKwall_glob(npro,9,nshl,nshl) 82 integer intp 83 84 85c 86c.... compute the nodes which lie on the boundary (hierarchic) 87c 88 call getbnodes(lnode) 89c 90c.... loop through the integration points 91c 92 if(lcsyst.eq.3.or.lcsyst.eq.4) then 93 ngaussb = nintb(lcsyst) 94 else 95 ngaussb = nintb(lcsyst) 96 endif 97 98 do intp = 1, ngaussb 99c 100c.... get the hierarchic shape functions at this int point 101c 102 call getshp(shpb, shglb, sgn, 103 & shape, shdrv) 104c 105c NOTE I DID NOT PASS THE lnode down. It is not needed 106c since the shape functions are zero on the boundary 107c 108c Note that xmudmi is not calculated at these quadrature 109c points so you give it a zero. This has implications. 110c the traction calculated by this approach will include 111c molecular stresses ONLY. This is why we will use the 112c consistent flux method to obtain the forces when doing 113c effective viscosity wall modeling. When doing slip velocity 114c this is not a problem since the traction is given from the 115c log law relation (not the viscosity). 116c 117 xmudmi=zero 118c 119c.... get necessary fluid properties (including eddy viscosity) 120c 121 call getdiff(dwl, yl, shape, xmudmi, xlb, rmu, rho) 122c 123c.... calculate the integraton variables 124c 125 call e3bvar (yl, acl, ul, 126 & shape, 127 & shdrv, xlb, 128 & lnode, WdetJb, 129 & bnorm, pres, 130 & u1, u2, u3, 131 & rmu, unm, 132 & tau1n, tau2n, tau3n, 133 & vdot, rlKwall, 134 & xKebe, rKwall_glob) 135 136c 137c.... -----------------> boundary conditions <------------------- 138c 139 do iel = 1, npro 140c 141c if we have a nonzero value then 142c calculate the fluxes through this surface 143c 144 iface = abs(iBCB(iel,2)) 145!MR CHANGE 146 if(iface>MAXSURF) then 147 write(*,*) 'iface>MAXSURF',iface,MAXSURF 148 write(*,*) 'Increase MAXSURF or decrease surfID in geom.spj' 149 stop !brutal but mpi_finalize will not work without broadcasting the information to other processors. 150 endif 151!MR CHANGE END 152 153 if (nsrflist(iface) .ne. 0 .and. ires.ne.2) then 154 flxID(1,iface) = flxID(1,iface) + WdetJb(iel)! measure area too 155 flxID(2,iface) = flxID(2,iface) - WdetJb(iel) * unm(iel) 156 flxID(3,iface) = flxID(3,iface) 157 & - ( tau1n(iel) - bnorm(iel,1)*pres(iel)) 158 & * WdetJb(iel) 159 flxID(4,iface) = flxID(4,iface) 160 & - ( tau2n(iel) - bnorm(iel,2)*pres(iel)) 161 & * WdetJb(iel) 162 flxID(5,iface) = flxID(5,iface) 163 & - ( tau3n(iel) - bnorm(iel,3)*pres(iel)) 164 & * WdetJb(iel) 165 166 endif 167c 168c 169c.... mass flux 170c 171 172 if (btest(iBCB(iel,1),0)) then 173 unm(iel) = zero 174 do n = 1, nshlb 175 nodlcl = lnode(n) 176 unm(iel) = unm(iel) 177 & + shape(iel,nodlcl) * BCB(iel,n,1) 178 enddo 179 endif 180c 181c.... pressure 182c 183 if (btest(iBCB(iel,1),1)) then 184 pres(iel) = zero 185 do n = 1, nshlb 186 nodlcl = lnode(n) 187 pres(iel) = pres(iel) 188 & + shape(iel,nodlcl) * BCB(iel,n,2) 189 enddo 190 endif 191c 192c.... viscous flux 193c 194 if (btest(iBCB(iel,1),2)) then 195 tau1n(iel) = zero 196 tau2n(iel) = zero 197 tau3n(iel) = zero 198 do n = 1, nshlb 199 nodlcl = lnode(n) 200 tau1n(iel) = tau1n(iel) 201 & + shape(iel,nodlcl)*BCB(iel,n,3) 202 tau2n(iel) = tau2n(iel) 203 & + shape(iel,nodlcl)*BCB(iel,n,4) 204 tau3n(iel) = tau3n(iel) 205 & + shape(iel,nodlcl)*BCB(iel,n,5) 206 enddo 207 endif 208!KJ CHANGE if(ideformwall.eq.1) then 209c 210c.... turbulence wall (as a way of checking for deformable wall stiffness) 211c 212 if (btest(iBCB(iel,1),4)) then 213 rlKwall(iel,:,:) = rlKwall(iel,:,:) / ngaussb ! divide by number of gauss points 214 pres(iel) = zero ! to avoid the gauss point loop 215 tau1n(iel) = zero ! and make the traction contribution 216 tau2n(iel) = zero ! zero 217 tau3n(iel) = zero 218 else 219 rlKwall(iel,:,:) = zero ! this is not a deformable element 220 vdot(iel,:) = zero 221 xKebe(iel,:,:,:) = zero 222 rKwall_glob(iel,:,:,:) = zero ! no stiffness: not a wall element 223 endif 224!KJ CHANGE endif 225 226 enddo ! end of bc loop 227c 228c$$$c.... if we are computing the bdry for the consistent 229c$$$c boundary forces, we must not include the surface elements 230c$$$c in the computataion (could be done MUCH more efficiently!)---> 231 !this 232 !comment should read as for the consistent flux calculation rather than boundary forces 233c$$$c 234 if (ires .eq. 2) then 235 do iel = 1, npro 236 if (nsrflist(iBCB(iel,2)) .ne. 0) then 237 unm(iel) = zero 238 tau1n(iel) = zero 239 tau2n(iel) = zero 240 tau3n(iel) = zero 241c pres(iel) = zero 242c 243c whatever is zeroed here will beome part of the post-processed surface 244c "traction force" 245!KJ CHANGE if(ideformwall.eq.1) then 246c 247c uncomment the next two lines to get all of the t vector coming from 248c Alberto's wall motion model. 249c vdot(iel,:)=zero 250c rlKwall(iel,:,:)=zero 251 252c 253c uncomment the next 8 lines to get only the tangential part 254c 255 vn=dot_product(vdot(iel,:),bnorm(iel,:)) 256 vdot(iel,:)=vn*bnorm(iel,:) 257 walln1=dot_product(rlkwall(iel,1,:),bnorm(iel,:)) 258 walln2=dot_product(rlkwall(iel,2,:),bnorm(iel,:)) 259 walln3=dot_product(rlkwall(iel,3,:),bnorm(iel,:)) 260 rlkwall(iel,1,:)=walln1*bnorm(iel,:) 261 rlkwall(iel,2,:)=walln2*bnorm(iel,:) 262 rlkwall(iel,3,:)=walln3*bnorm(iel,:) 263!KJ CHANGE endif 264 endif 265 enddo 266 endif 267c 268c.... assemble the contributions 269c 270 rNa(:,1) = -WdetJb * ( tau1n - bnorm(:,1) * pres - vdot(:,1)) 271 rNa(:,2) = -WdetJb * ( tau2n - bnorm(:,2) * pres - vdot(:,2)) 272 rNa(:,3) = -WdetJb * ( tau3n - bnorm(:,3) * pres - vdot(:,3)) 273 rNa(:,4) = WdetJb * unm 274c 275 if(iconvflow.eq.1) then ! conservative form was integrated 276 ! by parts and has a convective 277 ! boundary integral 278c 279c.... assemble the contributions 280c 281 rou=rho*unm 282 rNa(:,1) = rNa(:,1) + WdetJb * rou * u1 283 rNa(:,2) = rNa(:,2) + WdetJb * rou * u2 284 rNa(:,3) = rNa(:,3) + WdetJb * rou * u3 285 endif 286c 287c.... -------------------------> Residual <-------------------------- 288c 289c.... add the flux to the residual 290c 291 do n = 1, nshlb 292 nodlcl = lnode(n) 293 294 rl(:,nodlcl,1) = rl(:,nodlcl,1) - shape(:,nodlcl) * rNa(:,1) 295 rl(:,nodlcl,2) = rl(:,nodlcl,2) - shape(:,nodlcl) * rNa(:,2) 296 rl(:,nodlcl,3) = rl(:,nodlcl,3) - shape(:,nodlcl) * rNa(:,3) 297 rl(:,nodlcl,4) = rl(:,nodlcl,4) - shape(:,nodlcl) * rNa(:,4) 298 299 enddo 300 if(ideformwall.eq.1) then 301 rl(:,1,1) = rl(:,1,1) - rlKwall(:,1,1) 302 rl(:,1,2) = rl(:,1,2) - rlKwall(:,1,2) 303 rl(:,1,3) = rl(:,1,3) - rlKwall(:,1,3) 304 305 rl(:,2,1) = rl(:,2,1) - rlKwall(:,2,1) 306 rl(:,2,2) = rl(:,2,2) - rlKwall(:,2,2) 307 rl(:,2,3) = rl(:,2,3) - rlKwall(:,2,3) 308 309 rl(:,3,1) = rl(:,3,1) - rlKwall(:,3,1) 310 rl(:,3,2) = rl(:,3,2) - rlKwall(:,3,2) 311 rl(:,3,3) = rl(:,3,3) - rlKwall(:,3,3) 312 endif 313c 314c.... --------------------> Aerodynamic Forces <--------------------- 315c 316 if (((ires .ne. 2) .and. (iter .eq. nitr)) 317 & .and. (abs(itwmod).eq.1)) then 318c 319c.... compute the forces on the body 320c 321 where (nsrflist(iBCB(:,2)).eq.1) 322 tau1n = ( tau1n - bnorm(:,1)*pres) * WdetJb 323 tau2n = ( tau2n - bnorm(:,2)*pres) * WdetJb 324 tau3n = ( tau3n - bnorm(:,3)*pres) * WdetJb 325 elsewhere 326 tau1n = zero 327 tau2n = zero 328 tau3n = zero 329 endwhere 330c 331c Note that the sign has changed from the compressible code to 332c make it consistent with the way the bflux calculates the forces 333c Note also that Hflux has moved to e3btemp 334c 335 Force(1) = Force(1) - sum(tau1n) 336 Force(2) = Force(2) - sum(tau2n) 337 Force(3) = Force(3) - sum(tau3n) 338 339c 340 endif 341c 342c.... end of integration loop 343c 344 enddo 345 if(ideformwall.eq.1) then 346c 347c.... -----> Wall Stiffness and Mass matrices for implicit LHS <----------- 348c 349c.... Now we simply have to add the stiffness contribution in rKwall_glob to 350c.... the mass contribution already contained in xKebe 351 352c.... this line is going to destroy the mass matrix contribution 353 354 355c xKebe = zero 356 357 xKebe(:,:,:,:) = ( xKebe(:,:,:,:)*iwallmassfactor 358 & + rKwall_glob(:,:,:,:)*iwallstiffactor ) 359 360 endif 361c$$$ ttim(40) = ttim(40) + tmr() 362c 363c.... return 364c 365 return 366 end 367 368 369c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 370c********************************************************************* 371c********************************************************************* 372 373 374 subroutine e3bSclr (yl, iBCB, BCB, shpb, shglb, 375 & xlb, rl, sgn, dwl) 376 include "common.h" 377c 378 dimension yl(npro,nshl,ndof), iBCB(npro,ndiBCB), 379 & BCB(npro,nshlb,ndBCB), shpb(nshl,*), 380 & shglb(nsd,nshl,*), 381 & xlb(npro,nenl,nsd), 382 & rl(npro,nshl) 383c 384 real*8 WdetJb(npro), bnorm(npro,nsd) 385c 386 dimension lnode(27), sgn(npro,nshl), 387 & shape(npro,nshl), shdrv(npro,nsd,nshl), 388 & rNa(npro), flux(npro) 389 real*8 dwl(npro,nshl) 390 391c 392c.... compute the nodes which lie on the boundary (hierarchic) 393c 394 call getbnodes(lnode) 395c 396c.... loop through the integration points 397c 398 if(lcsyst.eq.3.or.lcsyst.eq.4) then 399 ngaussb = nintb(lcsyst) 400 else 401 ngaussb = nintb(lcsyst) 402 endif 403 do intp = 1, ngaussb 404c 405c.... get the hierarchic shape functions at this int point 406c 407 call getshp(shpb, shglb, sgn, 408 & shape, shdrv) 409c 410c.... calculate the integraton variables 411c 412 call e3bvarSclr (yl, shdrv, xlb, 413 & shape, WdetJb, bnorm, 414 & flux, dwl ) 415c 416c.... -----------------> boundary conditions <------------------- 417c 418 419c 420c.... heat or scalar flux 421c 422 if(isclr.eq.0) then 423 iwalljump=0 424 else 425 iwalljump=1 !turb wall between heat and scalar flux..jump over 426 endif 427 ib=4+isclr+iwalljump 428 ibb=6+isclr 429 do iel=1, npro 430c 431c if we have a nonzero value then 432c calculate the fluxes through this surface 433c 434 if (nsrflist(iBCB(iel,2)).ne.0 .and. ires.ne.2) then 435 iface = abs(iBCB(iel,2)) 436 flxID(ibb,iface) = flxID(ibb,iface) 437 & - WdetJb(iel) * flux(iel) 438 endif 439 440 if (btest(iBCB(iel,1),ib-1)) then 441 flux(iel) = zero 442 do n = 1, nshlb 443 nodlcl = lnode(n) 444 flux(iel) = flux(iel) 445 & + shape(iel,nodlcl) * BCB(iel,n,ibb) 446 enddo 447 endif 448 enddo 449c 450c.... assemble the contributions 451c 452 rNa(:) = -WdetJb * flux 453c 454c.... -------------------------> Residual <-------------------------- 455c 456c.... add the flux to the residual 457c 458 do n = 1, nshlb 459 nodlcl = lnode(n) 460 461 rl(:,nodlcl) = rl(:,nodlcl) - shape(:,nodlcl) * rNa(:) 462 enddo 463c 464c.... --------------------> Aerodynamic Forces <--------------------- 465c 466 if ((ires .ne. 2) .and. (iter .eq. nitr)) then 467c 468c.... compute the forces on the body 469c 470 if(isclr.eq.0) HFlux = sum(flux) 471c 472 endif 473c 474c.... end of integration loop 475c 476 enddo 477 478c 479c.... return 480c 481 return 482 end 483 484