From f503057793e4720ea22ced501f2d89ce82d5ca5b Mon Sep 17 00:00:00 2001 From: Timothy <75321887+timothy-nunn@users.noreply.github.com> Date: Mon, 2 Dec 2024 09:47:01 +0000 Subject: [PATCH] Convert fw.f90 to Python (#3387) * Convert friction and heat_transfer to Python * Convert fw_thermal_conductivity to Python --- CMakeLists.txt | 1 - process/blanket_library.py | 3 +- process/fw.py | 93 +++++++++++++++++- source/fortran/fw.f90 | 188 ------------------------------------- tests/unit/test_fw.py | 24 +++++ 5 files changed, 115 insertions(+), 194 deletions(-) delete mode 100644 source/fortran/fw.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 1af3fc15..a6b1de57 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -74,7 +74,6 @@ FILE(GLOB PROCESS_PYTHON_SRC_PATHS ${PROCESS_PYTHON_SRC_DIR}/*.py) LIST(APPEND PROCESS_SRCS numerics.f90 scan.f90 - fw.f90 hcpb.f90 pfcoil.f90 reinke_module.f90 diff --git a/process/blanket_library.py b/process/blanket_library.py index d3cb747d..c8a5004a 100644 --- a/process/blanket_library.py +++ b/process/blanket_library.py @@ -20,7 +20,6 @@ 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 @@ -2412,7 +2411,7 @@ def pressure_drop( # 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) + lamda = self.fw.friction(reyn) # Pressure drop coefficient diff --git a/process/fw.py b/process/fw.py index 084d5759..614d9955 100644 --- a/process/fw.py +++ b/process/fw.py @@ -4,8 +4,8 @@ constants, fwbs_variables, error_handling as eh, - fw_module, process_output as po, + error_handling, ) from process.utilities.f2py_string_patch import f2py_compatible_to_string from process.coolprop_interface import FluidProperties @@ -127,10 +127,10 @@ def fw_temp( eh.report_error(224) # Thermal conductivity of first wall material (W/m.K) - tkfw = fw_module.fw_thermal_conductivity(temp_k) + tkfw = self.fw_thermal_conductivity(temp_k) # Heat transfer coefficient (W/m2K) - hcoeff = fw_module.heat_transfer( + hcoeff = self.heat_transfer( masflx, ob_fluid_properties.density, afw, @@ -303,3 +303,90 @@ def fw_temp( ) return tpeakfw, cfmean, rhofmean, massrate + + def fw_thermal_conductivity(self, t): + """Calculates the thermal conductivity of the first wall + t : input real : property temperature (K) + Calculates the thermal conductivity of Eurofer (W/m/K). + """ + # Eurofer correlation, from "Fusion Demo Interim Structural Design Criteria - + # Appendix A Material Design Limit Data", F. Tavassoli, TW4-TTMS-005-D01, 2004 + # t in Kelvin + return ( + (5.4308 + 0.13565 * t - 0.00023862 * t * t + 1.3393e-7 * t * t * t) + * fwbs_variables.fw_th_conductivity + / 28.34 + ) + + def heat_transfer(self, masflx, rhof, radius, cf, viscf, kf): + """Calculate heat transfer coefficient using Gnielinski correlation + author: M Kovari, CCFE, Culham Science Centre + masflx : input real : coolant mass flux in a single channel (kg/m2/s) + rhof : input real : coolant density (average of inlet and outlet) (kg/m3) + radius : input real : coolant pipe radius (m) + cf : input real : coolant specific heat capacity (average of inlet and outlet) (J/K) + viscf : input real : coolant viscosity (average of inlet and outlet) (Pa.s) + kf : input real : thermal conductivity of coolant (average of inlet and outlet) (W/m.K) + Gnielinski correlation. Ignore the distinction between wall and + bulk temperatures. Valid for:3000 < Re < 5e6, 0.5 < Pr < 2000 + https://en.wikipedia.org/wiki/Nusselt_number#Gnielinski_correlation + """ + # Calculate pipe diameter (m) + diameter = 2 * radius + + # Calculate flow velocity (m/s) + velocity = masflx / rhof + + # Calculate Reynolds number + reynolds = rhof * velocity * diameter / viscf + + # Calculate Prandtl number + pr = cf * viscf / kf + + # Calculate Darcy friction factor, using Haaland equation + f = self.friction(reynolds) + + # Calculate the Nusselt number + nusselt = ( + (f / 8.0) + * (reynolds - 1000.0) + * pr + / (1 + 12.7 * np.sqrt(f / 8.0) * (pr**0.6667 - 1.0)) + ) + + # Calculate the heat transfer coefficient (W/m2K) + heat_transfer = nusselt * kf / (2.0 * radius) + + # Check that Reynolds number is in valid range for the Gnielinski correlation + if (reynolds <= 3000.0) or (reynolds > 5.0e6): + error_handling.fdiags[0] = reynolds + error_handling.report_error(225) + + # Check that Prandtl number is in valid range for the Gnielinski correlation + if (pr < 0.5) or (pr > 2000.0): + error_handling.fdiags[0] = pr + error_handling.report_error(226) + + # Check that the Darcy friction factor is in valid range for the Gnielinski correlation + if f <= 0.0: + error_handling.report_error(227) + + return heat_transfer + + def friction(self, reynolds): + """Calculate Darcy friction factor, using Haaland equation + author: M Kovari, CCFE, Culham Science Centre + reynolds : input real : Reynolds number + darcy_friction : output real : Darcy friction factor + Darcy friction factor, using Haaland equation, an approximation to the + implicit Colebrook-White equationGnielinski correlation. + https://en.wikipedia.org/wiki/Darcy_friction_factor_formulae#Haaland_equation + """ + + # Bracketed term in Haaland equation + bracket = ( + fwbs_variables.roughness / fwbs_variables.afw / 3.7 + ) ** 1.11 + 6.9 / reynolds + + # Calculate Darcy friction factor + return (1.8 * np.log10(bracket)) ** (-2) diff --git a/source/fortran/fw.f90 b/source/fortran/fw.f90 deleted file mode 100644 index ef5d2090..00000000 --- a/source/fortran/fw.f90 +++ /dev/null @@ -1,188 +0,0 @@ -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -module fw_module - !! Module containing first wall model - !! author: J Morris, CCFE, Culham Science Centre - !! N/A - !! This module contains the PROCESS first wall model - !! PROCESS Engineering paper (M. Kovari et al.) - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Modules to import -#ifndef dp - use, intrinsic :: iso_fortran_env, only: dp=>real64 -#endif - implicit none - - private - public :: friction, heat_transfer, fw_thermal_conductivity - -contains - - subroutine friction(reynolds, darcy_friction) - !! Calculate Darcy friction factor, using Haaland equation - !! author: M Kovari, CCFE, Culham Science Centre - !! reynolds : input real : Reynolds number - !! darcy_friction : output real : Darcy friction factor - !! Darcy friction factor, using Haaland equation, an approximation to the - !! implicit Colebrook–White equationGnielinski correlation. - !! https://en.wikipedia.org/wiki/Darcy_friction_factor_formulae#Haaland_equation - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use fwbs_variables, only: roughness, afw - - implicit none - - real(dp), intent(in) :: reynolds - real(dp), intent(out) :: darcy_friction - - ! Local variables ! - ! !!!!!!!!!!!!!!!!!! - - real(dp) :: bracket - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Bracketed term in Haaland equation - bracket = (roughness/afw/3.7)**1.11d0 + 6.9/reynolds - - ! Calculate Darcy friction factor - darcy_friction = (1.8d0 * log10(bracket))**(-2) - - end subroutine friction - - function heat_transfer(masflx, rhof, radius, cf, viscf, kf) - !! Calculate heat transfer coefficient using Gnielinski correlation - !! author: M Kovari, CCFE, Culham Science Centre - !! masflx : input real : coolant mass flux in a single channel (kg/m2/s) - !! rhof : input real : coolant density (average of inlet and outlet) (kg/m3) - !! radius : input real : coolant pipe radius (m) - !! cf : input real : coolant specific heat capacity (average of inlet and outlet) (J/K) - !! viscf : input real : coolant viscosity (average of inlet and outlet) (Pa.s) - !! kf : input real : thermal conductivity of coolant (average of inlet and outlet) (W/m.K) - !! Gnielinski correlation. Ignore the distinction between wall and - !! bulk temperatures. Valid for:3000 < Re < 5e6, 0.5 < Pr < 2000 - !! https://en.wikipedia.org/wiki/Nusselt_number#Gnielinski_correlation - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use error_handling, only: fdiags, report_error - - implicit none - - ! Arguments - - ! Function output: Heat transfer coefficient (W/m2K) - real(dp) :: heat_transfer - - ! Coolant mass flux in a single channel (kg/m2/s) - real(dp) :: masflx - - ! Coolant density (average of inlet and outlet) (kg/m3) - real(dp) :: rhof - - ! Coolant pipe radius (m) - real(dp) :: radius - - ! Coolant specific heat capacity (average of inlet and outlet) (J/K) - real(dp) :: cf - - ! Coolant viscosity (average of inlet and outlet) (Pa.s) - real(dp) :: viscf - - ! Thermal conductivity of coolant (average of inlet and outlet) (W/m.K) - real(dp) :: kf - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Local variables - - ! Calculate flow velocity (m/s) - real(dp) :: velocity - - ! Reynolds number - real(dp) :: reynolds - - ! Prandtl number - real(dp) :: pr - - ! Darcy friction factor - real(dp) :: f - - ! Nusselt number - real(dp) :: nusselt - - ! Pipe diameter (m) - real(dp) ::diameter - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Calculate pipe diameter (m) - diameter = 2*radius - - ! Calculate flow velocity (m/s) - velocity = masflx / rhof - - ! Calculate Reynolds number - reynolds = rhof * velocity * diameter / viscf - - ! Calculate Prandtl number - pr = cf * viscf / kf - - ! Calculate Darcy friction factor, using Haaland equation - call friction(reynolds, f) - - ! Calculate the Nusselt number - nusselt = (f/8.0d0)*(reynolds-1000.0d0)*pr / (1+12.7*sqrt(f/8.0d0)*(pr**0.6667-1.0d0)) - - ! Calculate the heat transfer coefficient (W/m2K) - heat_transfer = nusselt * kf / (2.0d0*radius) - - ! Check that Reynolds number is in valid range for the Gnielinski correlation - if ( ( reynolds <= 3000.0d0 ).or.( reynolds > 5.0d6 ) ) then - fdiags(1) = reynolds - call report_error(225) - end if - - ! Check that Prandtl number is in valid range for the Gnielinski correlation - if ( ( pr < 0.5d0 ).or.( pr > 2000.0d0) ) then - fdiags(1) = pr - call report_error(226) - end if - - ! Check that the Darcy friction factor is in valid range for the Gnielinski correlation - if ( f <= 0.0d0 ) call report_error(227) - - end function heat_transfer - - function fw_thermal_conductivity(t) - !! Calculates the thermal conductivity of the first wall - !! t : input real : property temperature (K) - !! Calculates the thermal conductivity of Eurofer (W/m/K). - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use fwbs_variables, only: fw_th_conductivity - - implicit none - - ! Function return value: thermal conductivity of first wall (W/m.K) - real(dp) :: fw_thermal_conductivity - - ! Arguments ! - ! !!!!!!!!!!!! - - real(dp), intent(in) :: t - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Eurofer correlation, from "Fusion Demo Interim Structural Design Criteria - - ! Appendix A Material Design Limit Data", F. Tavassoli, TW4-TTMS-005-D01, 2004 - ! t in Kelvin - fw_thermal_conductivity = (5.4308D0 + 0.13565D0*t - 0.00023862D0*t*t + & - 1.3393D-7*t*t*t)*fw_th_conductivity/28.34D0 - - end function fw_thermal_conductivity - -end module diff --git a/tests/unit/test_fw.py b/tests/unit/test_fw.py index 198cd3bc..36d8a6b0 100644 --- a/tests/unit/test_fw.py +++ b/tests/unit/test_fw.py @@ -161,3 +161,27 @@ def test_fw_temp(fwtempparam, monkeypatch, fw): assert rhofmean == pytest.approx(fwtempparam.expected_rhofmean, rel=1e-4) assert massrate == pytest.approx(fwtempparam.expected_massrate, rel=1e-4) + + +def test_friction(monkeypatch, fw): + monkeypatch.setattr(fwbs_variables, "afw", 0.1) + monkeypatch.setattr(fwbs_variables, "roughness", 1e-6) + + assert fw.friction(5500) == pytest.approx(0.0366668931278784) + + +def test_heat_transfer(fw): + assert fw.heat_transfer( + masflx=112.19853108876258, + rhof=8.8673250601290707, + radius=0.0060000000000000001, + cf=5184.9330299967578, + viscf=4.0416219836935569e-05, + kf=0.3211653052986152, + ) == pytest.approx(1929.2210192752084) + + +def test_fw_thermal_conductivity(monkeypatch, fw): + monkeypatch.setattr(fwbs_variables, "fw_th_conductivity", 28.9) + + assert fw.fw_thermal_conductivity(1900.0) == pytest.approx(326.70406785462256)