module blanket_library ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! This library contains routines that can be shared by the blanket modules used in PROCESS. !!! These include: !!! - component_volumes !!! - component_masses !!! - thermo_hydraulic_model !!! author: G Graham, CCFE, Culham Science Centre !! Acronyms for this module: !! !! BB Breeding Blanket !! FW First Wall !! BZ Breeder Zone !! MF/BSS Manifold/Back Supporting Structure !! LT Low Temperature !! HT High Temperature !! MMS Multi Module Segment !! SMS Single Modle Segment !! IB Inboard !! OB Outboard !! HCD Heating & Current Drive !! FCI Flow Channel Insert !!! Any changes within a subroutine or function code will have a comment explaining the change ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #ifndef dp use, intrinsic :: iso_fortran_env, only: dp=>real64 #endif implicit none real(dp) :: volshldi, volshldo !! Volume of inboard and outboard shield (m3) real(dp) :: volvvi, volvvo !! Volume of inboard and outboard Vacuum Vessel (m3) real(dp) :: hcryopf !! Clearance between uppermost PF coil and cryostat lid (m) real(dp) :: vfblkti, vfblkto !! Inboard/outboard void fraction of blanket real(dp) :: bldepti, bldepto !! Inboard/outboard blanket coolant channel length (radial direction) (m) real(dp) :: blwidti, blwidto !! Inboard/outboard blanket mid-plan toroidal circumference for segment (m) real(dp) :: bllengi, bllengo !! Inboard/outboard blanket segment poloidal length (m) real(dp) :: bzfllengi, bzfllengo !! Inboard/outboard primary blanket flow lengths (m) real(dp) :: bzfllengi_liq, bzfllengo_liq !! Inboard/outboard secondary blanket flow lengths (m) real(dp) :: pnucfwi, pnucfwo !! Inboard/outboard first wall nuclear heating (MW) real(dp) :: tpeakfwi, tpeakfwo !! Inboard/outboard first wall peak temperature (K) real(dp) :: mffwi, mffwo, mffw !! Inboard/outboard total mass flow rate to remove inboard FW power (kg/s) real(dp) :: npfwi, npfwo !! Inboard/utboard total number of pipes real(dp) :: mffwpi, mffwpo !! Inboard/outboard mass flow rate per coolant pipe (kg/s) real(dp) :: pnucblkti, pnucblkto !! Neutron power deposited inboard/outboard blanket blanket (MW) real(dp) :: mfblkti, mfblkto, mfblkt !! Inboard/outboard blanket mass flow rate for coolant (kg/s) real(dp):: mfblkti_liq, mfblkto_liq, mfblkt_liq !! Inboard/outboard blanket mass flow rate for liquid breeder (kg/s) real(dp) :: mftotal !! Total mass flow rate for coolant (kg/s) real(dp) :: npblkti, npblkto !! Inboard/outboard total num of pipes real(dp) :: mfblktpi, mfblktpo !! Inboard/outboard mass flow rate per coolant pipe (kg/s) real(dp) :: velblkti, velblkto !! Inboard/outboard coolant velocity in blanket (m/s) real(dp) :: htpmw_fwi, htpmw_fwo !! Inboard/outboard first wall pumping power (MW) real(dp) :: htpmw_blkti, htpmw_blkto !! Inboard/outboard blanket pumping power (MW) real(dp) :: htpmw_fw_blkti, htpmw_fw_blkto !! Inboard/outboard fw and blanket pumping power (MW) real(dp) :: hblnkt !! Blanket internal half-height (m) real(dp) :: hshld !! Shield internal half-height (m) real(dp) :: hvv !! Vacuum vessel internal half-height (m) integer :: icomponent !! Switch used to specify selected component: blanket=0, shield=1, vacuum vessel=2 contains subroutine init_blanket_library !! Initialise module variables implicit none hblnkt = 0.0D0 hshld = 0.0D0 hcryopf = 0.0D0 hvv = 0.0D0 volshldi = 0.0D0 volshldo = 0.0D0 volvvi = 0.0D0 volvvo = 0.0D0 bldepti = 0.0D0 bldepto = 0.0D0 blwidti = 0.0D0 blwidto = 0.0D0 bllengi = 0.0D0 bllengo = 0.0D0 bzfllengi = 0.0D0 bzfllengi_liq = 0.0D0 bzfllengo_liq = 0.0D0 bzfllengo = 0.0D0 pnucfwi = 0.0D0 pnucfwo = 0.0D0 tpeakfwi = 0.0D0 tpeakfwo = 0.0D0 mffwi = 0.0D0 mffwo = 0.0D0 mffw = 0.0D0 npfwi = 0.0D0 npfwo = 0.0D0 mffwpi = 0.0D0 mffwpo = 0.0D0 pnucblkti = 0.0D0 pnucblkto = 0.0D0 mfblkti = 0.0D0 mfblkto = 0.0D0 mfblkti_liq = 0.0D0 mfblkto_liq = 0.0D0 mfblkt_liq = 0.0D0 mfblkt = 0.0D0 mftotal = 0.0D0 npblkti = 0.0D0 npblkto = 0.0D0 mfblktpi = 0.0D0 mfblktpo = 0.0D0 velblkti = 0.0D0 velblkto = 0.0D0 htpmw_fwi = 0.0D0 htpmw_fwo = 0.0D0 htpmw_blkti = 0.0D0 htpmw_blkto = 0.0D0 vfblkti = 0.0D0 vfblkto = 0.0D0 end subroutine init_blanket_library ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! VOLUME CALCULATIONS !!! component_volumes: !!! - component_half_height !!! - dshaped_component !!! - elliptical_component !!! - apply_coverage_factors !!! - external_cryo_geometry ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine component_volumes !! Calculate the blanket, shield, vacuum vessel and cryostat volumes !! author: J. Morris, CCFE, Culham Science Centre !! Calculate the blanket, shield, vacuum vessel and cryostat volumes ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use fwbs_variables, only: fwbsshape use physics_variables, only: itart implicit none ! N.B. icomponent is a switch used to specify selected component: blanket=0, sheild=1, vacuum vessel=2 ! Replaced seperate subroutines for blnkt, shld and vv with fuction/subroutine with icomponent switch. ! Calculate half-height ! Blanket hblnkt = component_half_height(icomponent=0) ! Shield hshld = component_half_height(icomponent=1) ! Vacuum Vessel hvv = component_half_height(icomponent=2) ! D-shaped blanket and shield if ((itart == 1).or.(fwbsshape == 1)) then do icomponent = 0,2 call dshaped_component(icomponent) enddo ! Elliptical blanket and shield else do icomponent = 0,2 call elliptical_component(icomponent) enddo ! This will fail the hts_REBCO and 2D_scan regression tests, ! the number of VMCON iterations (nviter) is different. ! Seems to be because in the blanket calculations (icomponent=0): ! r2 = 1.3836567143743970 rather than old value of r2 = 1.3836567143743972, ! r3 = 3.7009701431231936 rather than r3 = 3.7009701431231923. end if ! Apply coverage factors to volumes and surface areas call apply_coverage_factors ! Calculate cryostat geometry call external_cryo_geometry end subroutine component_volumes ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function component_half_height(icomponent) !! Calculate the blanket, shield or vacuum vessel half-height !! Based on blanket_half_height, shield_half_height, vv_half_height !! original author: J. Morris, CCFE, Culham Science Centre !! author: G. Graham, CCFE, Culham Science Centre ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use build_variables, only: hmax, vgap_xpoint_divertor, vgap_vv_thermalshield, blnktth, shldtth, scrapli, scraplo, & fwith, fwoth, d_vv_bot, d_vv_top use physics_variables, only: rminor, kappa, idivrt use divertor_variables, only: divfix implicit none ! Input variables integer, intent(in) :: icomponent ! icomponent - blnkt=0, shld=1, vv=2 ! Local variables real(dp) :: hbot, htop ! Return variable real(dp) :: component_half_height ! Calculate component internal lower half-height (m) ! Blanket if (icomponent==0) hbot = rminor*kappa + vgap_xpoint_divertor + divfix - blnktth ! Sheild if (icomponent==1) hbot = rminor*kappa + vgap_xpoint_divertor + divfix ! Vacuum vessel if (icomponent==2) hbot = hmax - vgap_vv_thermalshield - d_vv_bot ! Calculate component internal upper half-height (m) ! If a double null machine then symmetric if (idivrt == 2) then htop = hbot else ! Blanket htop = rminor*kappa + 0.5D0*(scrapli+scraplo + fwith+fwoth) ! Shield if (icomponent==1) htop = htop + blnktth ! Vacuum Vessel if (icomponent==2) htop = htop + blnktth + shldtth end if ! Average of top and bottom (m) component_half_height = 0.5D0*(htop + hbot) end function component_half_height ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine dshaped_component(icomponent) !! Calculate component surface area and volume using dshaped scheme !! Based on dshaped_blanket, dshaped_shield, dshaped_vv !! original author: J. Morris, CCFE, Culham Science Centre !! author: G. Graham, CCFE, Culham Science Centre ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use maths_library, only: dshellarea, dshellvol use build_variables, only: rsldi, shldith, blnkith, fwith, scrapli, scraplo, & fwoth, blareaib, blareaob, blarea, blnkoth, blnktth, & shareaib, shareaob, sharea, shldoth, shldtth, & rsldo, d_vv_in, d_vv_out, d_vv_top, d_vv_bot use fwbs_variables, only: volblkti, volblkto, volblkt, volshld, vdewin use physics_variables, only: rminor implicit none ! Input variables integer, intent(in) :: icomponent ! icomponent - blnkt=0, shld=1, vv=2 ! Local variables real(dp) :: r1, r2 ! Calculate major radius to outer edge of inboard ... ! ... section (m) r1 = rsldi ! ... shield (m) if (icomponent==1) r1 = r1 + shldith ! ... blanket (m) if (icomponent==0) r1 = r1 + shldith + blnkith ! Horizontal distance between inside edges (m) ! i.e. outer radius of inboard part to inner radius of outboard part ! Blanket r2 = fwith + scrapli + 2.0D0*rminor + scraplo + fwoth ! Sheild if (icomponent==1) r2 = blnkith + r2 + blnkoth ! Vaccum Vessel if (icomponent==2) r2 = rsldo - r1 ! Calculate surface area, assuming 100% coverage if (icomponent==0) call dshellarea(r1, r2, hblnkt, blareaib, blareaob, blarea) if (icomponent==1) call dshellarea(r1, r2, hshld, shareaib, shareaob, sharea) ! Calculate volumes, assuming 100% coverage if (icomponent==0) call dshellvol(r1, r2, hblnkt, blnkith, blnkoth, blnktth, volblkti, volblkto, volblkt) if (icomponent==1) call dshellvol(r1, r2, hshld, shldith, shldoth, shldtth, volshldi, volshldo, volshld) if (icomponent==2) call dshellvol(r1, r2, hvv, d_vv_in, d_vv_out, & (d_vv_top+d_vv_bot)/2, volvvi, volvvo, vdewin) end subroutine dshaped_component ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine elliptical_component(icomponent) !! Calculate component surface area and volume using elliptical scheme !! Based on elliptical_blanket, elliptical_shield, elliptical_vv !! original author: J. Morris, CCFE, Culham Science Centre !! author: G. Graham, CCFE, Culham Science Centre ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use maths_library, only: eshellarea, eshellvol use build_variables, only: rsldi, shldith, blnkith, rsldo, shldoth, blnkoth, & blareaib, blareaob, blarea, blnktth, & shareaib, shareaob, sharea, shldtth, & d_vv_in, d_vv_out, d_vv_top, d_vv_bot use fwbs_variables, only: volblkti, volblkto, volblkt, volshld, vdewin use physics_variables, only: rmajor, rminor, triang implicit none ! Input variables integer, intent(in) :: icomponent ! icomponent - blnkt=0, shld=1, vv=2 ! Local variables real(dp) :: r1, r2, r3 ! Major radius to centre of inboard and outboard ellipses (m) ! (coincident in radius with top of plasma) r1 = rmajor - rminor*triang ! Calculate distance between r1 and outer edge of inboard ... ! ... section (m) r2 = r1 - rsldi ! ... shield (m) if (icomponent==1) r2 = r2 - shldith ! ... blanket (m) if (icomponent==0) r2 = r2 - shldith - blnkith ! Calculate distance between r1 and inner edge of outboard ... ! ... section (m) r3 = rsldo - r1 ! ... shield (m) if (icomponent==1) r3 = r3 - shldoth ! ... blanket (m) if (icomponent==0) r3 = r3 - shldoth - blnkoth ! Calculate surface area, assuming 100% coverage if (icomponent==0) call eshellarea(r1, r2, r3, hblnkt, blareaib, blareaob, blarea) if (icomponent==1) call eshellarea(r1, r2, r3, hshld, shareaib, shareaob, sharea) ! Calculate volumes, assuming 100% coverage if (icomponent==0) call eshellvol(r1, r2, r3, hblnkt, blnkith, blnkoth, blnktth, & volblkti, volblkto, volblkt) if (icomponent==1) call eshellvol(r1, r2, r3, hshld, shldith, shldoth, shldtth, & volshldi, volshldo, volshld) if (icomponent==2) call eshellvol(r1, r2, r3, hvv, d_vv_in, d_vv_out, & (d_vv_top+d_vv_bot)/2, volvvi, volvvo, vdewin) end subroutine elliptical_component ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine apply_coverage_factors !! Apply coverage factors to volumes !! author: J. Morris, CCFE, Culham Science Centre !! Apply coverage factors to volumes ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use build_variables, only: blareaob, blarea, blareaib, shareaib, shareaob, & sharea use fwbs_variables, only: fdiv, fhcd, volblkto, volblkti, volblkt, fvolsi, & fvolso, volshld, vdewin, fvoldw use physics_variables, only: idivrt implicit none ! Apply blanket coverage factors if (idivrt == 2) then ! double null configuration blareaob = blarea*(1.0D0-2.0D0*fdiv-fhcd) - blareaib else ! single null configuration blareaob = blarea*(1.0D0-fdiv-fhcd) - blareaib end if blarea = blareaib + blareaob volblkto = volblkt*(1.0D0-fdiv-fhcd) - volblkti volblkt = volblkti + volblkto ! Apply shield coverage factors shareaib = fvolsi*shareaib shareaob = fvolso*shareaob sharea = shareaib + shareaob volshldi = fvolsi*volshldi volshldo = fvolso*volshldo volshld = volshldi + volshldo ! Apply vacuum vessel coverage factor ! moved from dshaped_* and elliptical_* to keep coverage factor ! changes in the same location. vdewin = fvoldw*vdewin end subroutine apply_coverage_factors ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine external_cryo_geometry !! Calculate cryostat geometry !! author: J. Morris, CCFE, Culham Science Centre !! Calculate cryostat geometry ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use constants, only: pi use build_variables, only: clhsf, hmax, tfcth, ddwex use fwbs_variables, only: rdewex, rpf2dewar, zdewex, vdewex, vvmass, vdewin, & denstl, dewmkg use pfcoil_variables, only: rb, zh use buildings_variables, only: clh1 implicit none ! cryostat radius (m) ! ISSUE #508 Remove RFP option ! rb(i) = outer radius of PF coil i (tokamaks) rdewex = maxval(rb) + rpf2dewar ! Clearance between uppermost PF coil and cryostat lid (m). ! Scaling from ITER by M. Kovari hcryopf = clhsf * (2.0D0*rdewex)/28.440D0 ! Half-height of cryostat (m) ! ISSUE #508 Remove RFP option zdewex = maxval(zh) + hcryopf ! Vertical clearance between TF coil and cryostat (m) clh1 = zdewex - (hmax + tfcth) ! cryostat volume (m3) vdewex = ( (2.0D0*pi*rdewex) * 2.0D0*zdewex + (2.0D0*pi*rdewex**2) ) * ddwex ! Vacuum vessel mass (kg) vvmass = vdewin * denstl ! Sum of internal vacuum vessel and cryostat masses (kg) dewmkg = (vdewin + vdewex) * denstl end subroutine external_cryo_geometry ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! MASS CALCULATIONS !!! component_masses !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine component_masses !! Calculations for component masses !! author: J. Morris, CCFE, Culham Science Centre !! Calculations for component masses ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! variables used in ccfe hcpb but not kit hcpb: use constants, only: pi use build_variables, only: fwareaib, fwith, fwareaob, fwoth, fwarea use divertor_variables, only: divsur, divclfr, divplt, fdiva, divmas, divdens use physics_variables, only: rminor, rmajor, idivrt, sarea ! variables below used in kit hcpb but not ccfe hcpb: use build_variables, only: blnkith, blbuith, blbmith, blbpith, blnkoth, & blbuoth, blbmoth, blbpoth ! fwbs_variables: use fwbs_variables, only: volblkt, vfblkt, & !! CCFE and KIT HCPB -------- whtblbe, whtblss, denstl, whtblkt, & volshld, vfshld, coolmass, fwclfr, & !! CCFE HCPB only ------------------- breeder_f, breeder_multiplier, whtbltibe12, whtblli4sio4, wtblli2o, & vfcblkt, vfpblkt, whtshld, wpenshld, fwmass, fw_armour_vol, & fw_armour_thickness, fw_armour_mass, armour_fw_bl_mass, & volblkti, volblkto, iblnkith, fblhebmi, & !! KIT HCPB only --------------- fblhebpi,fblhebmo, fblhebpo, fblss, fblbe, & whtblbreed, densbreed, fblbreed, & iblanket, denw, vffwi, vffwo, volfw, & !! added -------------------- fblss_ccfe, fblli2sio4, fbltibe12 implicit none ! Only ccfe hcpb has local varibles i.e. coolvol ! Local variables ! Coolant volume (m3) real(dp) :: coolvol ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! CCFE HCPB modal calculates the coolant mass, ! have added an if staement using the iblanket switch for this. ! N.B. iblanket=1 for CCFE HCPB and iblanket=3 for the same with TBR using Shimwell. if ((iblanket==1).or.(iblanket==3)) then ! Start adding components of the coolant mass: ! Divertor coolant volume (m3) coolvol = divsur * divclfr * divplt ! Blanket coolant volume (m3) coolvol = coolvol + volblkt*vfblkt ! Shield coolant volume (m3) coolvol = coolvol + volshld*vfshld ! First wall coolant volume (m3) coolvol = coolvol + fwareaib*fwith*vffwi + fwareaob*fwoth*vffwo ! Mass of He coolant = volume * density at typical coolant temperatures and pressures (kg) coolmass = coolvol*1.517D0 ! Average first wall coolant fraction, only used by old routines in fispact.f90, safety.f90 fwclfr = (fwareaib*fwith*vffwi + fwareaob*fwoth*vffwo) / (fwarea*0.5D0*(fwith+fwoth)) endif ! CCFE HCPB calculates the mass of the divertor, blanket (including seprate masses for each material), ! shield, FW and FW armour. ! KIT HCPB calculates the mass of the blanket (including seprate masses for each material) ! and the void fraction for the blanket. ! N.B. iblanket=1 for CCFE HCPB and iblanket=3 for the same with TBR using Shimwell. if ((iblanket==1).or.(iblanket==3)) then ! Component masses ! Divertor mass (kg) divsur = fdiva * 2.0D0 * pi * rmajor * rminor if (idivrt == 2) divsur = divsur * 2.0D0 divmas = divsur * divdens * (1.0D0 - divclfr) * divplt ! Shield mass (kg) whtshld = volshld * denstl * (1.0D0 - vfshld) ! Penetration shield mass (set = internal shield) (kg) wpenshld = whtshld ! First wall volume (m^3) volfw = (fwareaib*fwith*(1.0D0-vffwi) + fwareaob*fwoth*(1.0D0-vffwo)) ! First wall mass, excluding armour (kg) fwmass = denstl * volfw ! First wall armour volume (m^3) fw_armour_vol = sarea*fw_armour_thickness ! First wall armour mass (kg) fw_armour_mass = fw_armour_vol*denw endif if ((iblanket==1).or.(iblanket==3)) then if (breeder_f < 1.0D-10) breeder_f = 1.0D-10 if (breeder_f > 1.0D0 ) breeder_f = 1.0D0 ! fbltibe12 = fblli2sio4 * (1 - breeder_f)/breeder_f ! New combined variable breeder_multiplier ! Lithium orthosilicate fraction: fblli2sio4 = breeder_f * breeder_multiplier ! Titanium beryllide fraction, and mass (kg): fbltibe12 = breeder_multiplier - fblli2sio4 whtbltibe12 = volblkt * fbltibe12 * 2260.0D0 ! Blanket Lithium orthosilicate mass (kg) ! Ref: www.rockwoodlithium.com... whtblli4sio4 = volblkt * fblli2sio4 * 2400.0D0 ! TODO sort this out so that costs model uses new variables. ! #327 For backwards compatibility, set the old blanket masses the same: whtblbe = whtbltibe12 wtblli2o = whtblli4sio4 ! Steel fraction by volume is the remainder: fblss_ccfe = 1.0D0 - fblli2sio4 - fbltibe12 - vfcblkt - vfpblkt ! Steel mass (kg) whtblss = volblkt * fblss_ccfe * denstl ! Total blanket mass (kg) whtblkt = whtbltibe12 + whtblli4sio4 + whtblss ! Total mass of first wall and blanket armour_fw_bl_mass = fw_armour_mass + fwmass + whtblkt endif if (iblanket==2) then ! Mass of steel in blanket (kg) if (iblnkith==1) then whtblss = denstl * ( volblkti/blnkith * ( blbuith * fblss + blbmith * (1.0D0-fblhebmi) + & blbpith * (1.0D0-fblhebpi) ) + volblkto/blnkoth * ( blbuoth * fblss + & blbmoth * (1.0D0-fblhebmo) + blbpoth * (1.0D0-fblhebpo) ) ) else whtblss = denstl * ( volblkto/blnkoth * ( blbuoth * fblss + & blbmoth * (1.0D0-fblhebmo) + blbpoth * (1.0D0-fblhebpo) ) ) end if ! Mass of beryllium in blanket (kg) if (iblnkith==1) then whtblbe = 1850.0D0 * fblbe * ( (volblkti * blbuith/blnkith) + & (volblkto * blbuoth/blnkoth) ) else whtblbe = 1850.0D0 * fblbe * (volblkto * blbuoth/blnkoth) end if ! Mass of breeder material in blanket (kg) if (iblnkith==1) then whtblbreed = densbreed * fblbreed * ( (volblkti * blbuith/blnkith) + & (volblkto * blbuoth/blnkoth) ) else whtblbreed = densbreed * fblbreed * (volblkto * blbuoth/blnkoth) end if ! Mass of blanket (kg) whtblkt = whtblss + whtblbe + whtblbreed ! Void fraction of blanket inboard portion if (iblnkith==1) then vfblkti = volblkti/volblkt * ( (blbuith/blnkith) * (1.0D0 - fblbe - fblbreed - fblss) & + (blbmith/blnkith) * fblhebmi + (blbpith/blnkith) * fblhebpi ) else vfblkti = 0.0D0 end if ! Void fraction of blanket outboard portion vfblkto = volblkto/volblkt * ( (blbuoth/blnkoth) * (1.0D0 - fblbe - fblbreed - fblss) & + (blbmoth/blnkoth) * fblhebmo + (blbpoth/blnkoth) * fblhebpo ) ! Void fraction of blanket vfblkt = vfblkti + vfblkto endif end subroutine component_masses ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! THERMOHYDRAUIC CALCULATIONS !!! thermo_hydraulic_model: !!! - blanket_mod_pol_height !!! - flow_velocity !!! - liquid_breeder_properties !!! - primary_coolant_properties !!! - hydraulic_diameter !!! - elbow_coeff !!! - pressure_drop !!! - liquid_breeder_pressure_drop_mhd !!! - deltap_tot !!! - pumppower !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine blanket_mod_pol_height !! Calculations for blanket module poloidal height !! author: J. Morris, CCFE, Culham Science Centre !! Calculations for blanket module poloidal height for D shaped and elliptical machines ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use constants, only: pi use build_variables, only: scrapli, scraplo use fwbs_variables, only: fwbsshape, nblktmodpi, fdiv, nblktmodpo, & iblanket use physics_variables, only: itart, rminor, idivrt, rmajor, triang implicit none ! Local variables ! Mid-plane distance from inboard to outboard side (m) real(dp) :: a ! Internal half-height of blanket (m) real(dp) :: b ! Calculate ellipse circumference using Ramanujan approximation (m) real(dp) :: ptor ! Major radius where half-ellipses 'meet' (m) real(dp) :: r1 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if ((itart == 1).or.(fwbsshape == 1)) then ! D-shaped machine ! Segment vertical inboard surface (m) bllengi = (2.0D0*hblnkt) / nblktmodpi ! Calculate perimeter of ellipse that defines the internal ! surface of the outboard first wall / blanket ! Mid-plane distance from inboard to outboard side (m) a = scrapli + 2.0D0*rminor + scraplo ! Internal half-height of blanket (m) b = hblnkt ! Calculate ellipse circumference using Ramanujan approximation (m) ptor = pi * ( 3.0D0*(a+b) - sqrt( (3.0D0*a + b)*(a + 3.0D0*b) ) ) ! Calculate blanket poloidal length and segment, subtracting divertor length (m) ! kit hcll version only had the single null option if (idivrt == 2) then ! Double null configuration bllengo = 0.5D0*ptor * (1.0D0 - 2.0D0*fdiv) / nblktmodpo else ! single null configuration bllengo = 0.5D0*ptor * (1.0D0 - fdiv) / nblktmodpo end if ! shape defined by two half-ellipses else ! Major radius where half-ellipses 'meet' (m) r1 = rmajor - rminor*triang ! Internal half-height of blanket (m) b = hblnkt ! Distance between r1 and nearest edge of inboard first wall / blanket (m) a = r1 - (rmajor - rminor - scrapli) ! Calculate ellipse circumference using Ramanujan approximation (m) ptor = pi * ( 3.0D0*(a+b) - sqrt( (3.0D0*a + b)*(a + 3.0D0*b) ) ) ! Calculate inboard blanket poloidal length and segment, subtracting divertor length (m) ! Assume divertor lies between the two ellipses, so fraction fdiv still applies ! kit hcll version only had the single null option if (idivrt == 2) then ! Double null configuration bllengi = 0.5D0*ptor * (1.0D0 - 2.0D0*fdiv) / nblktmodpi else ! single null configuration bllengi = 0.5D0*ptor * (1.0D0 - fdiv) / nblktmodpi end if ! Distance between r1 and inner edge of outboard first wall / blanket (m) a = rmajor + rminor + scraplo - r1 ! Calculate ellipse circumference using Ramanujan approximation (m) ptor = pi * ( 3.0D0*(a+b) - sqrt( (3.0D0*a + b)*(a + 3.0D0*b) ) ) ! kit hcll version only had the single null option ! Calculate outboard blanket poloidal length and segment, subtracting divertor length (m) if (idivrt == 2) then ! Double null configuration bllengo = 0.5D0*ptor * (1.0D0 - 2.0D0*fdiv) / nblktmodpo else ! single null configuration bllengo = 0.5D0*ptor * (1.0D0 - fdiv) / nblktmodpo end if end if end subroutine blanket_mod_pol_height ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function flow_velocity(i_channel_shape, mass_flow_rate, flow_density) !! Calculate the coolant flow velocity (m/s) for given pipe mass flow rate and pipe size/shape. !! N.B. Assumed that primary BB and FW coolants have same pipe radius (= afw). !! author: G. Graham, CCFE !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use fwbs_variables, only: afw, a_bz_liq, b_bz_liq use constants, only: pi implicit none !! Function return parameter !!!!! real(dp) :: flow_velocity !! Arguments !!!!!!!!!!!!!!!!!!!!! !! Swicth for circular or rectangular channel crossection. !! Shape depends on whether primary or secondary coolant. !! - =1 circle (primary) !! - =2 rectangle (secondary) integer, intent(in) :: i_channel_shape !! Coolant mass flow rate per pipe (kg/s) real(dp), intent(in) :: mass_flow_rate !! Coolant density real(dp), intent(in) :: flow_density !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! If primary coolant then circular channels assumed if (i_channel_shape==1) flow_velocity = mass_flow_rate / (flow_density*pi*afw*afw) !! If secondary coolant then rectangular channels assumed if (i_channel_shape==2) flow_velocity = mass_flow_rate / (flow_density * a_bz_liq * b_bz_liq) end function flow_velocity !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine liquid_breeder_properties(ip, ofile) !! Calculates the fluid properties of the Liquid Metal Breeder/Coolant in the Blanket BZ !! Uses middle value of input and output temperatures of Liquid Metal Breeder/Coolant !! Curently have PbLi but can expand with e.g., Lithium !! !! author: G Graham, CCFE !! !! References: !! !! [Mal1995] Malang and Mattas (1995), Comparison of lithium and the eutectic !! lead-lithium alloy, two candidate liquid metal breeder materials !! for self-cooled blankets, Fusion Engineering and Design 27, 399-406. !! !! [Mas2008] Mas de les Valles et al. (2008), Lead-lithium material database for !! nuclear fusion technology, Journal of Nuclear Materials, Vol. 376(6). !! !! [Mar2019] Martelli et al. (2019), Literature review of lead-lithium !! thermophysical properties, Fusion Engineering and Design, 138, 183-195. !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use fwbs_variables, only: inlet_temp_liq, outlet_temp_liq, a_bz_liq, b_bz_liq, den_liq, & specific_heat_liq, thermal_conductivity_liq, dynamic_viscosity_liq, electrical_conductivity_liq, & i_bb_liq, hartmann_liq, b_mag_blkt, iblnkith use physics_variables, only: bt, aspect, rmajor use build_variables, only: blnkith, blnkoth use error_handling, only: report_error implicit none !! Arguments !!!!!!!!!!!!!!!!!!!!! integer, intent(in) :: ip, ofile !! Local variables !!!!!!!!!!!!!!! !! Gas constant (J K-1 mol-1) real(dp) :: r !! mid temp of liquid metal (K) real(dp) :: mid_temp_liq !! Ratio of conductivity to dynamic viscosity real(dp) :: con_vsc_rat !! Array of valid temperature ranges for breeder propertites real(dp) :: t_ranges(5,2) !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Gas constant r = 8.314 !! Use mid temp if (inlet_temp_liq==outlet_temp_liq) then mid_temp_liq = outlet_temp_liq else mid_temp_liq = (inlet_temp_liq + outlet_temp_liq)*0.5 endif !! If the liquid metal is PbLi... if (i_bb_liq==0) then !! PbLi from [Mar2019] !! Constant pressure ~ 17 atmospheres ~ 1.7D6 Pa !! Li content is ~ 17% !! !! density kg m-3 T in Kelvin range = 508-880 K !! !! specific_heat J kg-1 K-1 T in Kelvin range = 508-880 K !! !! thermal_conductivity W m-1 K-1 T in Celcius range = 508-773 K !! !! dynamic_viscosity Pa s T in Celcius range = 508-873 K !! !! electrical_conductivity A V-1 m-1 T in Kelvin range = 600-800 K !! Caculate properties den_liq = 1.052D4*(1 - mid_temp_liq*1.13D-4) specific_heat_liq = 1.95D2 - mid_temp_liq*9.116D-3 thermal_conductivity_liq = 1.95 + (mid_temp_liq - 273.15)*1.96D-2 dynamic_viscosity_liq = 6.11D-3 -(2.257D-5 * (mid_temp_liq-273.15)) & + (3.766D-8 * (mid_temp_liq-273.15)**2) - (2.289D-11 * (mid_temp_liq-273.15)**3) t_ranges(1:4,1) = 508.0D0 t_ranges(1:4,2) = 880.0D0 electrical_conductivity_liq = 1.0D0/(1.03D-6 - (6.75D-11 * mid_temp_liq) + & (4.18D-13 * mid_temp_liq**2)) t_ranges(5,1) = 600.0D0 t_ranges(5,2) = 800.0D0 !! If the liquid metal is Li... else if (i_bb_liq==1) then !! Temporary - should be updated with information from Li reviews conducted at CCFE once completed !! Li Properties from [Mal1995] at 300 Celcius !! den_liq = 505 !! kg/m3 !! specific_heat_liq = 4260 !! J kg-1 K-1 !! thermal_conductivity_liq = 46 !! W m-1 K-1 !! dynamic_viscosity_liq = 1.0D-6 !! m2 s-1 !! electrical_conductivity_liq = 3.03D6 !! A V-1 m-1 !! New from 'Application of lithium in systems of fusion reactors. 1. Physical and chemical properties of lithium' !! Lyublinski et al., 2009, Plasma Devicec and Operations den_liq = 504.43D0 - (0.2729D0 * mid_temp_liq) - (8.0035D-5 * mid_temp_liq**2) + (3.799D-8 * mid_temp_liq**3) specific_heat_liq = 31.227 + (0.205D6 * mid_temp_liq**(-2)) - (5.265D-3 * mid_temp_liq) + (2.628D6 * mid_temp_liq**(-2)) !! thermal_conductivity_liq also in paper dynamic_viscosity_liq = exp(-4.16D0 - (0.64D0 * log(mid_temp_liq)) + (262.1/mid_temp_liq)) electrical_conductivity_liq = (0.9249D9 * mid_temp_liq) + 2.3167D6 - (0.7131D3 * mid_temp_liq) endif !! Magnetic feild strength in T for Hartmann calculation !! IB if (iblnkith==1) b_mag_blkt(1) = bt * rmajor/(rmajor - (rmajor/aspect) - (blnkith/2)) !! We do not use this if there is no IB blanket, but will use edge as fill value if (iblnkith==0) b_mag_blkt(1) = bt * rmajor/(rmajor - (rmajor/aspect)) !! OB b_mag_blkt(2) = bt * rmajor/(rmajor + (rmajor/aspect) + (blnkoth/2)) !! Calculate Hartmann number con_vsc_rat = electrical_conductivity_liq/dynamic_viscosity_liq !! Use toroidal width of the rectangular cooling channel as characteristic length scale hartmann_liq = b_mag_blkt * a_bz_liq/2.0D0 * sqrt(con_vsc_rat) !! Error for temperature range of breeder property realtions if (i_bb_liq == 0) then if ((any(t_ranges(:,1) > mid_temp_liq)).or.(any(t_ranges(:,2) < mid_temp_liq))) call report_error(280) endif !! Output !!!!!!!!!!!!!!!!!!!!!!!! if (ip == 0) return call write_output_liquid_breeder_properties contains subroutine write_output_liquid_breeder_properties use process_output, only: oheadr, osubhd, ovarrf, ovarre, & ocmmnt, ovarin, ovarst use fwbs_variables, only: i_bb_liq, den_liq, dynamic_viscosity_liq, & electrical_conductivity_liq, icooldual implicit none !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call oheadr(ofile,'Blanket : Liquid Breeder Properties') if (icooldual == 1) call ocmmnt(ofile, 'Single coolant: liquid metal circulted for tritium extraction.') if (icooldual == 2) call ocmmnt(ofile, 'Dual coolant: self-cooled liquid metal breeder.') if (i_bb_liq == 0) call ocmmnt(ofile, 'Blanket breeder type (i_bb_liq=0), PbLi (~ 17% Li)') if (i_bb_liq == 1) call ocmmnt(ofile, 'Blanket breeder type (i_bb_liq=1), Li') call ovarrf(ofile, 'Density (kg m-3)', '(den_liq)', den_liq, 'OP ') call ovarrf(ofile, 'Viscosity (Pa s)', '(dynamic_viscosity_liq)', dynamic_viscosity_liq, 'OP ') call ovarrf(ofile, 'Electrical Conductivity (A V-1 m-1)', '(electrical_conductivity_liq)', electrical_conductivity_liq, 'OP ') call ovarrf(ofile, 'Hartmann Number IB', '(hartmann_liq)', hartmann_liq(1), 'OP ') call ovarrf(ofile, 'Hartmann Number OB', '(hartmann_liq)', hartmann_liq(2), 'OP ') call ovarre(ofile, 'Inlet Temperature (Celcius)', '(inlet_temp_liq)', inlet_temp_liq, 'OP ') call ovarre(ofile, 'Outlet Temperature (Celcius)', '(outlet_temp_liq)', outlet_temp_liq, 'OP ') if (i_bb_liq == 0) then if ((any(t_ranges(:,1) > mid_temp_liq)).or.(any(t_ranges(:,2) < mid_temp_liq))) then call ocmmnt(ofile, 'Outside temperature limit for one or more liquid metal breeder properties.') call ovarrf(ofile, 'Liquid metal temperature (K)', '(mid_temp_liq)', mid_temp_liq, 'OP ') call ocmmnt(ofile, 'Density: Max T = 880 K, Min T = 508 K') call ocmmnt(ofile, 'Specific heat: Max T = 880 K, Min T = 508 K') call ocmmnt(ofile, 'Thermal conductivity: Max T = 880 K, Min T = 508 K') call ocmmnt(ofile, 'Dynamic viscosity : Max T = 880 K, Min T = 508 K') call ocmmnt(ofile, 'Electrical conductivity: Max T = 800 K, Min T = 600 K') endif endif end subroutine write_output_liquid_breeder_properties end subroutine liquid_breeder_properties !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function pressure_drop(ip, ofile, i_ps, num_90, num_180, l_pipe, den, vsc, vv, label) !! Pressure drops are calculated for a pipe with a number of 90 !! and 180 degree bends. The pressure drop due to frictional forces along !! the total straight length of the pipe is calculated, then the pressure !! drop due to the bends is calculated. The total pressure drop is the sum !! of all contributions. !! !! original author: P. J. Knight, CCFE !! moved from previous version of pumppower function by: G Graham, CCFE !! !! N.B Darcy-Weisbach Equation (straight pipe): !! !! kstrght = lambda * L/D !! !! pressure drop = kstrght * (rho*V^2)/2 !! !! lambda - Darcy friction factor, L - pipe length, D - hydraulic diameter, !! rho - fluid density, V - fluid flow average velocity !! !! This function also calculates pressure drop equations for elbow bends, !! with modified coefficients. !! !! N.B. Darcy friction factor is estimated from the Haaland approximation. !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use fwbs_variables, only: afw, a_bz_liq, b_bz_liq use fw_module, only: friction implicit none !! Function return parameter !!!!! !! Pressure drop along the pipe (Pa) real(dp) :: pressure_drop !! Arguments !!!!!!!!!!!!!!!!!!!!! integer, intent(in) :: ip, ofile !! Swicth for primary or secondary coolant !! - =1 primary !! - =2 secondary integer, intent(in) :: i_ps !! Number of 90 and 180 degree bends in pipe integer, intent(in) :: num_90, num_180 !! Total flow length along pipe (m) real(dp), intent(in) :: l_pipe !! Coolant density (kg/m3) real(dp), intent(in) :: den !! Coolant viscosity (Pa s) real(dp), intent(in) :: vsc !! Coolant flow velocity (m/s) real(dp), intent(in) :: vv !! Description of this calculation character(len=*), intent(in) :: label !! Local variables !!!!!!!!!!!!!!! !! Hydraulic diameter of coolant flow channels (m) real(dp) :: dh !! Reynolds number real(dp) :: reyn !! Darcy friction factor real(dp) :: lambda !! Pressure drops for straight setions, 90 bends and 180 bends (Pa) real(dp) :: pdropstraight, pdrop90, pdrop180 !! Elbow radius for 90 and 180 deg brends (m) real(dp) :: elbow_radius !! Pressure drop coefficients real(dp) :: kelbwn, kelbwt, kstrght !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Calculate hydraulic dimater for round or retancular pipe (m) dh = hydraulic_diameter(i_ps) !! Reynolds number reyn = den * vv * dh / vsc !! Calculate Darcy friction factor !! N.B. friction function Uses Haaland approx. which assumes a filled circular pipe. !! Use dh which allows us to do fluid calculations for non-cicular tubes !! (dh is estimate appropriate for fully developed flow). call friction(reyn,lambda) !! Pressure drop coefficient !!!!! !! Straight section kstrght = lambda * l_pipe/dh !! In preveious version of pumppower: !! - elbow radius assumed = 0.018m for 90 degree elbow, from WCLL !! - elbow radius assumed half that of 90 deg case for 180 deg elbow !! Intialised value for afw is 0.006m, so elbow radius = 3 * afw, !! aka 1.5 * pipe diameter, which seems to be engineering standard for !! a steel pipe long-radius elbow (short-radius elbow = 2 * afw). !! If primary coolant... if (i_ps==1) then elbow_radius = 3 * afw !! If secondary coolant... else !! See DCLL elbow_radius = b_bz_liq endif !! 90 degree elbow pressure drop coefficient kelbwn = elbow_coeff(elbow_radius, 90.0D0, lambda, dh) !! 180 degree elbow pressure drop coefficient kelbwt = elbow_coeff(elbow_radius/2, 180.0D0, lambda, dh) !! Total (Pa) pdropstraight = kstrght * 0.5D0*den*vv*vv pdrop90 = num_90*kelbwn * 0.5D0*den*vv*vv pdrop180 = num_180*kelbwt * 0.5D0*den*vv*vv pressure_drop = pdropstraight + pdrop90 + pdrop180 if (ip==0) return call write_output_pressure_drop contains subroutine write_output_pressure_drop use process_output, only: oheadr, osubhd, ovarrf, ovarre, & ocmmnt, ovarin, ovarst, oblnkl use global_variables, only: verbose implicit none !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call osubhd(ofile, 'Pressure drop (friction) for ' // label) call ovarre(ofile, 'Reynolds number', '(reyn)', reyn, 'OP ') call ovarre(ofile, 'Darcy friction factor', '(lambda)', lambda, 'OP ') call ovarre(ofile, 'Pressure drop (Pa)', '(pressure_drop)', pressure_drop, 'OP ') call ocmmnt(ofile, 'This is the sum of the following:') call ovarre(ofile, ' Straight sections (Pa)', '(pdropstraight)', & pdropstraight, 'OP ') call ovarre(ofile, ' 90 degree bends (Pa)', '(pdrop90)', pdrop90, 'OP ') call ovarre(ofile, ' 180 degree bends (Pa)', '(pdrop180)', pdrop180, 'OP ') call ocmmnt(ofile, 'Additional information is printed when verbose = 1') if (verbose==1) then call oblnkl(ofile) call ovarre(ofile, 'Straight section pressure drop coefficient', & '(kstrght)', kstrght, 'OP ') call ovarre(ofile, '90 degree elbow coefficient', '(kelbwn)', kelbwn, 'OP ') call ovarre(ofile, '180 degree elbow coefficient coefficient', '(kelbwt)', kelbwt, 'OP ') end if end subroutine write_output_pressure_drop end function pressure_drop !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function hydraulic_diameter(i_channel_shape) !! Caculate the hydraulic diameter (m) for a given coolant pipe size/shape. !! author: G. Graham !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use fwbs_variables, only: afw, a_bz_liq, b_bz_liq implicit none !! Function return parameter !!!!! real(dp) :: hydraulic_diameter !! Arguments !!!!!!!!!!!!!!!!!!!!! !! Swicth for circular or rectangular channel crossection !! Shape depends on whether primary or secondary coolant !! - =1 circle (primary) !! - =2 rectangle (secondary) integer, intent(in) :: i_channel_shape !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! If primary coolant then circular channels assumed if (i_channel_shape==1) hydraulic_diameter = 2.0D0*afw !! If secondary coolant then rectangular channels assumed if (i_channel_shape==2) hydraulic_diameter = 2*a_bz_liq*b_bz_liq/(a_bz_liq+b_bz_liq) end function hydraulic_diameter !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function elbow_coeff(r_elbow, ang_elbow, lambda, dh) !! Function calculates elbow bends coefficients for pressure drop !! calculations. !! !! author: G. Graham, CCFE !! !! References: !! !! [Ide1969] Idel'Cik, I. E. (1969), Memento des pertes de charge, !! Collection de la Direction des Etudes et Recherches d'Electricité de France. !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use constants, only: pi implicit none !! Function return parameter !!!!! !! Elbow coefficient for pressure drop calculation real(dp) :: elbow_coeff !! Arguments !!!!!!!!!!!!!!!!!!!!! !! Pipe elbow radius (m) real(dp), intent(in) :: r_elbow !! Pipe elbow angle (degrees) real(dp), intent(in) :: ang_elbow !! Darcy Friction Factor real(dp), intent(in) :: lambda !! Hydraulic Diameter (m) real(dp), intent(in) :: dh !! Local variables !!!!!!!!!!!!!!! !! A and B from singularity coefficient equation real(dp) :: a, b !! Singularity and friction coefficients real(dp) :: xifn, xift, ximn, ximt !! Ratio of elbow radius to Darcy friction factor real(dp) :: r_ratio !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Pressure loss for bends calculated using equivilent lengths. !! coeff = lambda * L/D --> get equiv. of L if (ang_elbow==90) then a = 1.0D0 else if (ang_elbow<70) then a = 0.9D0 * sin(ang_elbow * pi/180.0D0) else if (ang_elbow>100) then a = 0.7D0 + (0.35D0 * sin((ang_elbow/90.0D0) * (pi/180.0D0))) else write(*,*) 'No formula for 70 <= elbow angle(deg) <= 100, only 90 deg option available in this range.' stop 1 endif r_ratio = r_elbow/dh if (r_ratio > 1) b = 0.21D0/r_ratio**0.5 if (r_ratio < 1) b = 0.21D0/r_ratio**2.5 if (r_ratio == 1) b = 0.21D0 !! Singularity ximt = a * b !! Friction xift = 0.0175D0 * lambda * (r_elbow/dh) * ang_elbow !! Elbow Coefficient elbow_coeff = ximt + xift end function elbow_coeff !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function liquid_breeder_pressure_drop_mhd(ip, ofile, vel, vsc, conduct_liq, l_channel, num_pol, label) !! Calculates the pressure drop in a liquid metal flow channel due to MHD effects. The total pressure !! drop is the sum of contributions. This is only used for secondary coolant/breeder so rectangular flow !! channels are assumed. !! !! author: G Graham, CCFE !! !! References: !! !! [Miy1986] Miyazaki et al. (1986), Magneto-Hydro-Dynamic Pressure Drop of Lithium !! Flow in Rectangular Ducts, Fusion Technology, 10:3P2A, 830-836, DOI: 10.13182/FST10-830 !! !! [Mal1995] Malang and Mattas (1995), Comparison of lithium and the eutectic !! lead-lithium alloy, two candidate liquid metal breeder materials !! for self-cooled blankets, Fusion Engineering and Design 27, 399-406 !! !! [Iba2013] Ibano et al (2013), Nutronics and pumping power analysis on the !! Tokamak reactor for the fusion-biomass hybrid concept, !! Fusion Engineering and Design, 88 !! !! [Sho2018] Shoki et al (2018), MHD pressure drop measurement of PbLi flow !! in double-bended pipe, Fusion Engineering and Design, 136, 17-23 !! !! [Klu2019] Kluber et al. (2019), Numerical simulations of 3D magnetohydrodynamic !! flows in dual-coolant lead lithium blankets, Fusion Engineering and Design, !! 146, 684-687 !! !! [Sua2021] MHD effects in geometrical sigularities on high velocity breeding !! blanket designs. Part II, ENR-PRD.BB-T007-D002, EFDA_D_2PDT9U. !! Also, see asssociated paper: Suarez et al. (2021), On the use of CFD !! to obtain head loss coefficients in hydraulic systems and it's appliaction !! to liquid metal flows in nuclear fusion reactor blankets, Plasma. Phys. !! Control fusion, 63, 124002 !! !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use fwbs_variables, only: ifci, a_bz_liq, b_bz_liq, hartmann_liq, b_mag_blkt, & bz_channel_conduct_liq, th_wall_secondary use constraint_variables, only: maxradwallload use physics_variables, only: btot implicit none !! Function return parameter !!!!! real(dp) :: liquid_breeder_pressure_drop_mhd !! Arguments !!!!!!!!!!!!!!!!!!!!! integer, intent(in) :: ip, ofile !! Liquid metal coolant/breeder flow velocity (m/s) real(dp), intent(in) :: vel !! Liquid metal visosity real(dp), intent(in) :: vsc !! Liquid metal conductivity real(dp), intent(in) :: conduct_liq !! Length long poloidal sections of channel real(dp), intent(in) :: l_channel !! Number long poloidal sections of channel integer, intent(in) :: num_pol !! Description of this calculation character(len=*), intent(in) :: label !! Local Variables !!!!!!!!!!!!!!! !! Half-widths of channel (m) real(dp) :: half_wth_a, half_wth_b !! MHD pressure drop for single channel real(dp) :: mhd_pressure_drop !! Magnetic field strenght (T) real(dp) :: b_mag !! Internal resistance of fluid real(dp) :: r_i !! Resistance of side walls real(dp) :: r_w !! Resistance of electrode walls real(dp) :: r_e !! Wall-to-fluid conductance ratio real(dp) :: big_c !! Term calculated for [Miy1986] pressure drop real(dp) :: kp !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !! Magnetic feild strength in IB or OB blanket if (label=='Inboard blanket breeder liquid') b_mag = b_mag_blkt(1) !! IB if (label=='Outboard blanket breeder liquid') b_mag = b_mag_blkt(2) !! OB !! Half-widths !! N.B. a_bz_liq (width in the toroidal direction) is in B direction half_wth_a = a_bz_liq * 0.5 half_wth_b = b_bz_liq *0.5 !! If have thin conducting walls... if (ifci/=1) then !! Caculate resistances of fluid and walls r_i = half_wth_b/(conduct_liq * half_wth_a) r_w = half_wth_b/(bz_channel_conduct_liq * th_wall_secondary) r_e = half_wth_a/(bz_channel_conduct_liq * th_wall_secondary) big_c = r_i/r_w !! Calculate pressure drop for conducting wall [Miy1986] kp = big_c/(1 + half_wth_a/(3 * half_wth_b) + big_c) mhd_pressure_drop = kp * conduct_liq * vel * (b_mag**2) * l_channel !! If have perfcetly insulating FCIs... else !! Calculate pressure drop for (perfectly) insulating FCI [Mal1995] mhd_pressure_drop = vel * b_mag * l_channel * sqrt(conduct_liq * vsc / half_wth_a) endif !! Total (Pa) liquid_breeder_pressure_drop_mhd = num_pol * mhd_pressure_drop if (ip==0) return call write_output_liquid_breeder_pressure_drop_mhd contains subroutine write_output_liquid_breeder_pressure_drop_mhd use process_output, only: oheadr, osubhd, ovarrf, ovarre, & ocmmnt, ovarin, ovarst implicit none !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! call osubhd(ofile, 'Liquid metal breeder/coolant MHD pressure drop for ' // label) if (ifci==0) then call ocmmnt(ofile, 'Flow channels have thin conducting walls (ifci==0)') call ovarre(ofile, 'Wall conductance (A V-1 m-1)', '(bz_channel_conduct_liq)', bz_channel_conduct_liq, 'OP ') else if (ifci==2) then call ocmmnt(ofile, 'Flow Channel Inserts (FCIs) used (ifci==2)') call ovarre(ofile, 'FCI conductance (A V-1 m-1)', '(bz_channel_conduct_liq)', bz_channel_conduct_liq, 'OP ') else call ocmmnt(ofile, 'Flow Channel Inserts - assumed perfect insulator (ifci==1)') endif call ovarre(ofile, 'Length of long poloidal secion of channel (m)', '(l_channel)', l_channel, 'OP ') call ovarin(ofile, 'Number of long poloidal secions of channel', '(num_pol)', num_pol, 'OP ') call ovarre(ofile, 'MHD pressure drop (Pa)', '(liquid_breeder_pressure_drop_mhd)', liquid_breeder_pressure_drop_mhd, 'OP ') end subroutine write_output_liquid_breeder_pressure_drop_mhd end function liquid_breeder_pressure_drop_mhd !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end module blanket_library