From 8c38f883ea786160c936426a649a482c0f83cfe9 Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Fri, 8 Nov 2024 09:46:00 +0000 Subject: [PATCH] Convert pressure_drop to Python --- process/blanket_library.py | 176 ++++++++++++++++++++++++++--- source/fortran/blanket_library.f90 | 148 ------------------------ tests/unit/test_blanket_library.py | 73 ++++++++++++ 3 files changed, 234 insertions(+), 163 deletions(-) diff --git a/process/blanket_library.py b/process/blanket_library.py index 9acd7475..903977c1 100644 --- a/process/blanket_library.py +++ b/process/blanket_library.py @@ -20,6 +20,7 @@ pfcoil_variables, buildings_variables, error_handling, + fw_module, ) from process.utilities.f2py_string_patch import f2py_compatible_to_string from process.coolprop_interface import FluidProperties @@ -493,7 +494,7 @@ def primary_coolant_properties(self, output: bool): "Calculated using mid temp(s) of system (or systems if use different collant types).", ) - # FW (or FW/BB)!!!!!!!!!!!!!!!!!! + # FW (or FW/BB) if fwbs_variables.ipump == 1: po.osubhd(self.outfile, "First Wall :") @@ -543,7 +544,7 @@ def primary_coolant_properties(self, output: bool): "OP ", ) - # BB !!!!!!!!!!!!!!!!!!!!!!!!!!!! + # BB if fwbs_variables.ipump == 1: po.osubhd(self.outfile, "Breeding Blanket :") @@ -1270,11 +1271,11 @@ def liquid_breeder_properties(self, output: bool = False): # 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 + # 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 @@ -1923,7 +1924,7 @@ def thermo_hydraulic_model(self, output: bool): self.outfile, "Summary of first wall and blanket thermohydraulics" ) - # FW !!!!!!!!!!!!!!!!!!!!!!!!!!!! + # FW po.osubhd(self.outfile, "First wall: ") po.ovarst( @@ -1990,7 +1991,7 @@ def thermo_hydraulic_model(self, output: bool): "OP ", ) - # BB !!!!!!!!!!!!!!!!!!!!!!!!!!!! + # BB po.osubhd(self.outfile, "Breeding Blanket (primary): ") po.ovarin( @@ -2036,7 +2037,7 @@ def thermo_hydraulic_model(self, output: bool): "OP ", ) - # BB Liquid Metal Breeder !!!!!!! + # BB Liquid Metal Breeder ! if fwbs_variables.icooldual > 0: po.osubhd(self.outfile, "Breeding Blanket (breeder): ") @@ -2082,7 +2083,7 @@ def thermo_hydraulic_model(self, output: bool): "OP ", ) - # Pumping Power !!!!!!!!!!!!!!!!!!!!!!!! + # Pumping Power po.osubhd(self.outfile, "Mechanical pumping power: ") if fwbs_variables.ipump == 1: @@ -2169,9 +2170,7 @@ def deltap_tot( :param no180: Number of 180 degree bends in pipe """ # Friction - for all coolants - frict_drop = blanket_library.pressure_drop( - int(output), - self.outfile, + frict_drop = self.pressure_drop( icoolpump, no90, no180, @@ -2180,6 +2179,7 @@ def deltap_tot( coolant_dynamic_viscosity, flow_velocity, label, + output=output, ) if icoolpump == 2: @@ -2216,6 +2216,152 @@ def deltap_tot( return deltap_tot + def pressure_drop( + self, + i_ps: int, + num_90: float, + num_180: float, + l_pipe: float, + den: float, + vsc: float, + vv: float, + label: str, + output: bool = False, + ): + """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 + + :param i_ps: switch for primary or secondary coolant + :param num_90: number of 90 degree bends in the pipe + :param num_180: number of 180 degree bends in the pipe + :param l_pipe: total flow length along pipe (m) + :param den: coolant density (kg/m3) + :param vsc: coolant viscosity (Pa s) + :param vv: coolant flow velocity (m/s) + :param label: component name + :param output: boolean of whether to write data to output file + + 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. + """ + + # Calculate hydraulic dimater for round or retancular pipe (m) + dh = blanket_library.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). + lamda = fw_module.friction(reyn) + + # Pressure drop coefficient !!!!! + + # Straight section + kstrght = lamda * 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: + elbow_radius = 3 * fwbs_variables.afw + # If secondary coolant... + else: + # See DCLL + elbow_radius = fwbs_variables.b_bz_liq + + # 90 degree elbow pressure drop coefficient + kelbwn = blanket_library.elbow_coeff(elbow_radius, 90.0, lamda, dh) + + # 180 degree elbow pressure drop coefficient + kelbwt = blanket_library.elbow_coeff(elbow_radius / 2, 180.0, lamda, dh) + + # Total (Pa) + pdropstraight = kstrght * 0.5 * den * vv * vv + pdrop90 = num_90 * kelbwn * 0.5 * den * vv * vv + pdrop180 = num_180 * kelbwt * 0.5 * den * vv * vv + + pressure_drop = pdropstraight + pdrop90 + pdrop180 + + if output: + po.osubhd(self.outfile, f"Pressure drop (friction) for {label}") + po.ovarre(self.outfile, "Reynolds number", "(reyn)", reyn, "OP ") + po.ovarre(self.outfile, "Darcy friction factor", "(lambda)", lamda, "OP ") + po.ovarre( + self.outfile, + "Pressure drop (Pa)", + "(pressure_drop)", + pressure_drop, + "OP ", + ) + po.ocmmnt(self.outfile, "This is the sum of the following:") + po.ovarre( + self.outfile, + " Straight sections (Pa)", + "(pdropstraight)", + pdropstraight, + "OP ", + ) + po.ovarre( + self.outfile, + " 90 degree bends (Pa)", + "(pdrop90)", + pdrop90, + "OP ", + ) + po.ovarre( + self.outfile, + " 180 degree bends (Pa)", + "(pdrop180)", + pdrop180, + "OP ", + ) + + # TN: always write verbose stuff, it has no harm + po.ovarre( + self.outfile, + "Straight section pressure drop coefficient", + "(kstrght)", + kstrght, + "OP ", + ) + po.ovarre( + self.outfile, "90 degree elbow coefficient", "(kelbwn)", kelbwn, "OP " + ) + po.ovarre( + self.outfile, + "180 degree elbow coefficient coefficient", + "(kelbwt)", + kelbwt, + "OP ", + ) + + return pressure_drop + def pumppower( self, output, @@ -2247,7 +2393,7 @@ def pumppower( :param primary_coolant_switch: Switch for FW/blanket coolant, (1=He or 2=H2O) if icoolpump=1 :param coolant_density: Density of coolant or liquid breeder """ - # Pumping power !!!!!!!!!!!!!!!!! + # Pumping power ! # Outlet pressure is 'pressure' # Inlet pressure (Pa) diff --git a/source/fortran/blanket_library.f90 b/source/fortran/blanket_library.f90 index 9a742767..b1e69e1f 100644 --- a/source/fortran/blanket_library.f90 +++ b/source/fortran/blanket_library.f90 @@ -197,154 +197,6 @@ end subroutine init_blanket_library !!! - pumppower !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine pressure_drop(ip, ofile, i_ps, num_90, num_180, l_pipe, den, vsc, vv, label, pressure_drop_out) - - !! 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, roughness - use fw_module, only: friction - use process_output, only: oheadr, osubhd, ovarrf, ovarre, & - ocmmnt, ovarin, ovarst - - implicit none - - !! Function return parameter !!!!! - - !! Pressure drop along the pipe (Pa) - real(dp), intent(out) :: pressure_drop_out - - !! 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) - call hydraulic_diameter(i_ps, dh) - - !! 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 - call elbow_coeff(elbow_radius, 90.0D0, lambda, dh, kelbwn) - - !! 180 degree elbow pressure drop coefficient - call elbow_coeff(elbow_radius/2, 180.0D0, lambda, dh, kelbwt) - - !! 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_out = pdropstraight + pdrop90 + pdrop180 - - if (ip==0) return - - 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_out, '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') - - end subroutine pressure_drop - !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine hydraulic_diameter(i_channel_shape, hydraulic_diameter_out) diff --git a/tests/unit/test_blanket_library.py b/tests/unit/test_blanket_library.py index db9f5100..9347d320 100644 --- a/tests/unit/test_blanket_library.py +++ b/tests/unit/test_blanket_library.py @@ -1661,3 +1661,76 @@ def test_liquid_breeder_properties( assert fwbs_variables.b_mag_blkt == pytest.approx( liquidbreederpropertiesparam.expected_b_mag_blkt ) + + +class PressureDropParam(NamedTuple): + afw: Any = None + a_bz_liq: Any = None + b_bz_liq: Any = None + roughness: Any = None + ip: Any = None + ofile: Any = None + i_ps: Any = None + num_90: Any = None + num_180: Any = None + l_pipe: Any = None + den: Any = None + vsc: Any = None + vv: Any = None + label: Any = None + expected_pressure_drop_out: Any = None + + +@pytest.mark.parametrize( + "pressuredropparam", + ( + PressureDropParam( + afw=0.0060000000000000001, + a_bz_liq=0.20000000000000001, + b_bz_liq=0.20000000000000001, + roughness=9.9999999999999995e-07, + ip=0, + ofile=11, + i_ps=1, + num_90=2, + num_180=0, + l_pipe=4, + den=10.405276820718059, + vsc=3.604452999475736e-05, + vv=32.753134225223164, + label="Inboard first wall", + expected_pressure_drop_out=36214.869527556766, + ), + ), +) +def test_pressure_drop(pressuredropparam, monkeypatch, blanket_library_fixture): + """ + Automatically generated Regression Unit Test for pressure_drop. + + This test was generated using data from blanket_files/large_tokamak_primary_pumping2.IN.DAT. + + :param pressuredropparam: the data used to mock and assert in this test. + :type pressuredropparam: pressuredropparam + + :param monkeypatch: pytest fixture used to mock module/class variables + :type monkeypatch: _pytest.monkeypatch.monkeypatch + """ + monkeypatch.setattr(fwbs_variables, "afw", pressuredropparam.afw) + monkeypatch.setattr(fwbs_variables, "a_bz_liq", pressuredropparam.a_bz_liq) + monkeypatch.setattr(fwbs_variables, "b_bz_liq", pressuredropparam.b_bz_liq) + monkeypatch.setattr(fwbs_variables, "roughness", pressuredropparam.roughness) + + pressure_drop_out = blanket_library_fixture.pressure_drop( + i_ps=pressuredropparam.i_ps, + num_90=pressuredropparam.num_90, + num_180=pressuredropparam.num_180, + l_pipe=pressuredropparam.l_pipe, + den=pressuredropparam.den, + vsc=pressuredropparam.vsc, + vv=pressuredropparam.vv, + label=pressuredropparam.label, + ) + + assert pressure_drop_out == pytest.approx( + pressuredropparam.expected_pressure_drop_out + )