1 module blowerControl 2 logical :: BC_enable 3 integer :: nBlower 4 real*8 :: BC_dt 5 real*8 :: BC_istep0 6 real*8 :: BC_t 7 8 type :: blowerData 9 integer :: surfID, !surface ID to use for blower inlet 10 & nmap, !number of inlet nodes 11 & mode !constant = 0, trapezoid = 1, sinusoid = 2 12 integer, allocatable :: map(:) !mapping to the inlet nodes 13 real*8, allocatable :: vscale(:) !scalling map used to create an inlet boundary layer 14 logical :: updateBC, !true if the BC dynamically updated 15 & enable !true to enable blower on with the particular surface ID 16 17 !The wave form is expected to look something like: 18 ! |--------- (a) ---------| 19 ! ________ _______ _ vmax 20 ! / \ / | 21 ! / \ / | 22 !____/ \__________/ |_ vmin 23 ! |-|--------|-| 24 ! (b) (c) (d) 25 26 !Temporal parameters 27 real*8 :: t_cycle !period of the cycle (a) 28 real*8 :: t_riseTime, t_fallTime !time of rising edge (b), falling edge (d) 29 real*8 :: t_fullOn !duration pulse is on full (c) 30 real*8 :: t0 !t = 0 reference for wave form 31 32 !State values 33 real*8 :: vmax, vmin !maximum, minimum blower velocity 34 real*8 :: T !uniform temperature at blower inlet 35 real*8 :: nu !uniform eddy viscosity at blower inlet 36 real*8 :: deltaBL !Boundary layer thickness 37 real*8 :: deltaBLsclr !Boundary layer thickness for scalar 38 39 end type blowerData 40 41 type(blowerData), allocatable :: blower(:) 42 43 contains 44 subroutine BC_setVars(nBlowers, blowerMode, 45 & surfID, enable, 46 & t_cycle, t_fullOn, 47 & t_riseTime, t_fallTime, 48 & vmax, vmin, 49 & T, nu, 50 & deltaBL, deltaBLsclr ) 51 52 integer :: nBlowers 53 integer :: blowerMode(nBlowers), 54 & surfID( nBlowers), enable( nBlowers) 55 real*8 :: t_cycle( nBlowers), t_fullOn( nBlowers), 56 & t_riseTime(nBlowers), t_fallTime( nBlowers), 57 & vmax( nBlowers), vmin( nBlowers), 58 & T( nBlowers), nu( nblowers), 59 & deltaBL( nBlowers), deltaBLsclr(nBlowers) 60 61 nBlower = nBlowers 62 63 if(.not. allocated(blower)) then 64 allocate(blower(nBlower)) 65 endif 66 67 do i = 1, nBlower 68 blower(i)%mode = blowerMode(i) 69 blower(i)%surfID = surfID(i) 70 blower(i)%enable = (enable(i) == 1) 71 blower(i)%t_cycle = t_cycle(i) 72 blower(i)%t_fullOn = t_fullOn(i) 73 blower(i)%t_riseTime = t_riseTime(i) 74 blower(i)%t_fallTime = t_fallTime(i) 75 blower(i)%vmax = vmax(i) 76 blower(i)%vmin = vmin(i) 77 blower(i)%T = T(i) 78 blower(i)%nu = nu(i) 79 blower(i)%deltaBL = deltaBL(i) 80 blower(i)%deltaBLsclr = deltaBLsclr(i) 81 enddo 82 83 end subroutine 84 85 86 subroutine BC_iter(BC) 87 !Updates BC with the correct velocities for the blower. 88 !Note: itrBC still needs to be called after blowerIter to update 89 ! values in yold. 90 91! use blowerControl 92 use wallData 93 include "common.h" 94 95 integer :: i, iBlower, nNodeMap 96 real*8 :: vmagref, vmag 97 real*8 :: t 98 real*8 :: BC(nshg, ndofBC) 99 100 !local variables to make the code shorter 101 real*8 :: t0, t1, t2, t3, t4, vmax, vmin 102 103 if(.not. BC_enable) return 104 105 BC_t = BC_dt*(istep - BC_istep0) !Update time 106 107 do iBlower = 1, nBlower 108 if(.not. blower(iBlower)%enable .or. 109 & .not. blower(iBlower)%updateBC ) cycle 110 111 !------------------------ 112 !Setup usefull variables 113 !------------------------ 114 vmax = blower(iBlower)%vmax 115 vmin = blower(iBlower)%vmin 116 117 if(blower(iBlower)%mode == 1) then !Trapezoid 118 !The wave form is expected to look something like: 119 ! ________ _______ _ vmax 120 ! / \ / | 121 ! / \ / | 122 !____/ \__________/ |_ vmin 123 ! | | | | | 124 ! t0 t1 t2 t3 t4 125 126 t0 = 0.0d00 127 t1 = t0 + blower(iBlower)%t_riseTime 128 t2 = t1 + blower(iBlower)%t_fullOn 129 t3 = t2 + blower(iBlower)%t_fallTime 130 t4 = blower(iBlower)%t_cycle 131 132 !Calculate the velocity magnitude 133 !t = mod(BC_dt*istep + blower(iBlower)%t0, t4) 134 t = mod(BC_t - blower(iBlower)%t0, t4) 135 136 if(t0 <= t .and. t < t1) then !rising edge 137 vmagref = (vmax - vmin)*(t - t0)/(t1 - t0) + vmin 138 elseif(t1 <= t .and. t < t2) then !high 139 vmagref = vmax 140 elseif(t2 <= t .and. t < t3) then !falling edge 141 vmagref = (vmax - vmin)*(t - t3)/(t2 - t3) + vmin 142 elseif(t3 <= t .and. t < t4) then !low 143 vmagref = vmin 144 endif 145 elseif(blower(iBlower)%mode == 2) then !sinusoid 146 ! _.--._ _.--._ |- vmax 147 ! _-' '-_ _-' ' | 148 !__..-' '-..__..-' |- vmin 149 ! | 150 ! | |------ t_cycle ------| 151 ! t0 152 153 !t is number of waves since t0 154 t = (BC_t - blower(iBlower)%t0)/blower(iBlower)%t_cycle 155 vmagref = 0.5*(vmax - vmin)*sin(2*pi*t) + 0.5*(vmax + vmin) 156 elseif(blower(iBlower)%mode == 0) then !Constant 157 vmagref = max(vmax, vmin) 158 endif 159 160 !---------- 161 !update BC 162 !---------- 163 do i = 1,blower(iBlower)%nmap 164 imapped = blower(iBlower)%map(i) 165 vmag = vmagref*blower(iBlower)%vscale(i) !scale the reference 166 !to create a BL 167 BC(imapped, 3) = -vmag*wnorm(imapped, 1) !x velocity 168 BC(imapped, 4) = -vmag*wnorm(imapped, 2) !y velocity 169 BC(imapped, 5) = -vmag*wnorm(imapped, 3) !z velocity 170 enddo !end loop over blower surface nodes 171 enddo 172 173 end subroutine 174 175 176 subroutine BC_Finalize() 177 !Deallocates allocatable variables in blower and then deallocates blower. 178 179 integer :: iBlower 180 181 if(BC_enable) then 182 do iBlower = 1,nBlower 183 deallocate(blower(iBlower)%map) 184 deallocate(blower(iBlower)%vscale) 185 enddo 186 endif 187 188 deallocate(blower) !note that blower gets allocated regardless 189 !of whether BC_enable is true 190 end subroutine 191 192 193 subroutine BC_ReadPhase(tstep) 194 include "common.h" 195 include "mpif.h" 196 197 integer :: tstep 198 integer :: i, nBlowerIn, surfIDin 199 character(len=128) :: fname 200 logical :: existFname 201 real*8 :: phi 202 203 if(.not. BC_enable) return !nothing needs to be written 204 205 BC_t = BC_dt*(istep - istep0) !Update time 206 if(myrank == master) then 207 !Even if blowerPhase.%i.dat exists, it may not have the same 208 !number of blowers; zero all t0 in blower before proceeding. 209 do i = 1,nBlower 210 blower(i)%t0 = 0.0d00 211 enddo 212 213 !Assemble the file name and test if it exists. 214 write(fname, "('blowerPhase.', I0, '.dat')") tstep 215 inquire(file=fname, exist=existFname) 216 217 if(existFname) then 218 open(unit=1003, file=fname, status='old') 219 220 read(1003, *) nBlowerIn 221 222 !In the future, the values should probably be matched by 223 !surface ID, but I don't have time for that now. 224 do i = 1,nBlowerIn 225 read(1003, *) surfIDin, phi 226 blower(i)%t0 = -phi*blower(i)%t_cycle 227 !Note the negative sign; t0 is the reference starting 228 !time. For positive phase, the cycle started sometime in 229 !the past. 230 enddo 231 232 close(1003) 233 endif !existFname 234 endif !myrank == master 235 236 !At this point, t0 on rank 0 should be correct; 237 !broadcast the contents of the file 238 if(numpe > 1) then 239 do i = 1,nBlower 240 call MPI_BCAST(blower(i)%t0, 1, MPI_DOUBLE_PRECISION, 241 & master, MPI_COMM_WORLD, ierr) 242 enddo 243 endif 244 end subroutine 245 246 247 subroutine BC_writePhase(tstep) 248 include "common.h" 249 250 integer :: tstep !time step to use in the file name 251 252 character(len=128) :: fname 253 integer :: i 254 real*8 :: phi 255 256 if(.not. BC_enable) return !nothing needs to be written 257 258 BC_t = BC_dt*(istep - BC_istep0) !Update time 259 260 if(myrank == master) then 261 write(fname, "('blowerPhase.', I0, '.dat')") tstep 262 open(unit=1003, file=fname, status='unknown') 263 264 write(1003, "(i9)") nBlower 265 do i = 1,nBlower 266 !Calculate the present possition in the blowing cycle 267 if(blower(i)%updateBC) then 268 phi = mod((BC_t - blower(i)%t0)/blower(i)%t_cycle,1.0d0) 269 else 270 phi = 0.0d00 271 endif 272 273 write(1003, "(i9, F18.14)") blower(i)%surfID, phi 274 enddo 275 276 close(1003) 277 endif 278 end subroutine 279 280! subroutine BC_setNBlower(nBlowers) 281! integer :: nBlowers 282! nBlower = nBlowers 283! end subroutine 284! 285! subroutine BC_setSurfID(surfID) 286! integer :: surfID(nBlower) 287! do i = 1, nBlower 288! blower(i)%surfID = surfID(i) 289! enddo 290! end subroutine 291! 292! subroutine BC_setEnable(enable) 293! integer :: enable(nBlower) 294! do i = 1, nBlower 295! blower(i)%enable = enable(i) 296! enddo 297! end subroutine 298! 299! subroutine BC_setSurfID(surfID) 300! real*8 :: surfID(nBlower) 301! do i = 1, nBlower 302! blower(i)%surfID = surfID(i) 303! enddo 304! end subroutine 305! 306! subroutine BC_setSurfID(surfID) 307! integer :: surfID(nBlower) 308! do i = 1, nBlower 309! blower(i)%surfID = surfID(i) 310! enddo 311! end subroutine 312 313 314 end module 315 316 subroutine BC_init(dt, tstep, BC) 317 use blowerControl 318 use wallData !wnorm, 319 use turbSA !d2wall 320 include "common.h" 321 322 real*8 :: BC(nshg, ndofBC) 323 324 integer :: i, iBlower, imapped, tstep 325 real*8 :: deltaBLinv, vmag, dt 326 327 !local variables to shorten the code a bit 328 integer :: nNodeMap 329 real*8 :: vmax !T, nu 330 331 !variables that should be set in solver.inp: 332 !nBlower + surfID, t_cycle, t_riseTime, t_fallTime, t_fullOn, vmax, vmin, T, nu, deltaBL, enable 333 334 !Test if any of the blowers are enabled. If not, then set 335 !BC_enable to false and exit. 336 BC_enable = .false. 337 do iBlower = 1, nBlower 338 if(blower(iBlower)%enable) then 339 BC_enable = .true. 340 exit 341 endif 342 enddo 343 if(.not. BC_enable) return 344 345 BC_istep0 = istep 346 BC_dt = dt 347 348 !check to see if a blowerControl.dat has been written with phase data 349 call BC_ReadPhase(tstep) !t0 is updated in the call. 350 351 do iBlower = 1, nBlower 352 !------------------------------------ 353 !Nodal Mapping and Referrence Arrays 354 !------------------------------------ 355 356 !create the nodal mapping 357 allocate(blower(iBlower)%map(nshg)) 358 call sfID2np(blower(iBlower)%surfID, !analog to sfID2np, except 359 & blower(iBlower)%nmap, !that map is allocated in 360 & blower(iBlower)%map) !the subroutine call. 361! call asfID2np(blower(iBlower)%surfID, !analog to sfID2np, except 362! & blower(iBlower)%nmap, !that map is allocated in 363! & blower(iBlower)%map) !the subroutine call. 364 nNodeMap = blower(iBlower)%nmap 365 366 !Calculate the boundary layer scalling 367 allocate(blower(iBlower)%vscale(nNodeMap)) 368 369 if(blower(iBlower)%deltaBL > 0.0) then !use the provided BL 370 deltaBLinv = 1/blower(iBlower)%deltaBL 371 else !use a default of essentially zero thickness 372 deltaBLinv = 1e16 373 endif 374 375 do i = 1, nNodeMap 376 !loop over each surface node and calculate vscale with a 377 !linear ramp. Add more options in the future. 378 imapped = blower(iBlower)%map(i) 379 blower(iBlower)%vscale(i) = 380 & min(1.0d0, d2wall(imapped)*deltaBLinv) 381 enddo 382 383 !------------------------------------- 384 !Apply persistent boundary conditions 385 !------------------------------------- 386 387 !temperature 388 if(blower(iBlower)%T > 0.0) then 389 do i = 1, nNodeMap 390 imapped = blower(iBlower)%map(i) 391 BC(imapped, 2) = blower(iBlower)%T 392 enddo 393 endif 394 395 !eddy viscosity 396 if(blower(iBlower)%nu >= 0.0) then 397 if(blower(iBlower)%deltaBLsclr > 0.0) then !use the provided BL 398 deltaBLinv = 1/blower(iBlower)%deltaBLsclr 399 else !use a default of essentially zero thickness 400 deltaBLinv = 1e16 401 endif 402 403 do i = 1, nNodeMap 404 imapped = blower(iBlower)%map(i) 405 BC(imapped, 7) = blower(iBlower)%nu * 406 & min(1.0d0, d2wall(imapped)*deltaBLinv) 407 enddo 408 endif 409 410 !velocity 411 blower(iBlower)%updateBC = .true. !update velocity later with a call to BC_iter. 412 enddo !loop over blower surfaces 413 414 call BC_iter(BC) 415 416 do iBlower = 1, nBlower 417 if(blower(iBlower)%t_cycle >= 1.0 .or. !Unreasonable cycle 418 & blower(iBlower)%t_cycle <= 0.0 .or. !periods are provided, 419 & blower(iBlower)%mode == 0) then !so assume const. 420 blower(iBlower)%updateBC = .false. 421 else 422 blower(iBlower)%updateBC = .true. 423 endif 424 enddo 425 426 end subroutine 427 428 429