diff --git a/src/Finch_Inputs.hpp b/src/Finch_Inputs.hpp index e30895a..92971a8 100644 --- a/src/Finch_Inputs.hpp +++ b/src/Finch_Inputs.hpp @@ -91,8 +91,11 @@ struct Properties { double density; double specific_heat; - double thermal_conductivity; - double thermal_diffusivity; + double thermal_conductivity_solid_0; + double thermal_conductivity_solid_T; + double thermal_conductivity_liquid_0; + double thermal_conductivity_liquid_T; + double vaporization_temperature; double latent_heat; double solidus; double liquidus; @@ -211,8 +214,12 @@ class Inputs Info << "Properties:" << std::endl; Info << " Density: " << properties.density << std::endl; Info << " Specific Heat: " << properties.specific_heat << std::endl; - Info << " Thermal Conductivity: " << properties.thermal_conductivity - << std::endl; + Info << " Thermal Conductivity, Solid: " + << properties.thermal_conductivity_solid_0 << " + " + << properties.thermal_conductivity_solid_T << " * T" << std::endl; + Info << " Thermal Conductivity, Liquid: " + << properties.thermal_conductivity_liquid_0 << " + " + << properties.thermal_conductivity_liquid_T << " * T" << std::endl; Info << " Latent Heat: " << properties.latent_heat << std::endl; Info << " Solidus: " << properties.solidus << std::endl; Info << " Liquidus: " << properties.liquidus << std::endl; @@ -270,13 +277,16 @@ class Inputs write(); - // create auxiliary properties - properties.thermal_diffusivity = - ( properties.thermal_conductivity ) / + // create auxiliary properties - estimate time step based on a maximum + // temperature and the parsed thermal conducitivity values + const double est_max_thermal_diffusivity = + ( properties.thermal_conductivity_liquid_0 + + properties.thermal_conductivity_liquid_T * + properties.vaporization_temperature ) / ( properties.density * properties.specific_heat ); time.time_step = ( time.Co * space.cell_size * space.cell_size ) / - ( properties.thermal_diffusivity ); + ( est_max_thermal_diffusivity ); Info << "Calculated time step: " << time.time_step << std::endl; @@ -335,8 +345,37 @@ class Inputs // Read properties components properties.density = db["properties"]["density"]; properties.specific_heat = db["properties"]["specific_heat"]; - properties.thermal_conductivity = - db["properties"]["thermal_conductivity"]; + // Thermal conductvity as a function of temperature in the form a + b * + // Temperature for solid and liquid. If only one value is given, fixed + // thermal conductivity for both phases are used + if ( db["properties"].contains( "thermal_conductivity_solid" ) && + db["properties"].contains( "thermal_conductivity_liquid" ) ) + { + properties.thermal_conductivity_solid_0 = + db["properties"]["thermal_conductivity_solid"][0]; + properties.thermal_conductivity_solid_T = + db["properties"]["thermal_conductivity_solid"][1]; + properties.thermal_conductivity_liquid_0 = + db["properties"]["thermal_conductivity_liquid"][0]; + properties.thermal_conductivity_liquid_T = + db["properties"]["thermal_conductivity_liquid"][1]; + // Also expects a max temperature (vaporization temperature) for + // thermal conductivity calculations + properties.vaporization_temperature = + db["properties"]["vaporization_temperature"]; + } + else + { + properties.thermal_conductivity_solid_0 = + db["properties"]["thermal_conductivity"]; + properties.thermal_conductivity_liquid_0 = + db["properties"]["thermal_conductivity"]; + properties.thermal_conductivity_solid_T = 0.0; + properties.thermal_conductivity_liquid_T = 0.0; + // Default value for a max temperature (unused since + // thermal_conductivity_T = 0) + properties.vaporization_temperature = 5000.0; + } properties.latent_heat = db["properties"]["latent_heat"]; properties.solidus = db["properties"]["solidus"]; properties.liquidus = db["properties"]["liquidus"]; diff --git a/src/Finch_Solver.hpp b/src/Finch_Solver.hpp index 5653edb..d0457fe 100644 --- a/src/Finch_Solver.hpp +++ b/src/Finch_Solver.hpp @@ -44,9 +44,15 @@ class Solver double dt_; double solidus_; double liquidus_; + double inv_freezing_range_; double rho_cp_; double rho_Lf_by_dT_; - double k_by_dx2_; + double inv_dx2_; + double k_solid_0_; + double k_solid_T_; + double k_liquid_0_; + double k_liquid_T_; + double temp_max_; // heat source parameters double power_; @@ -73,11 +79,20 @@ class Solver liquidus_ = db.properties.liquidus; + inv_freezing_range_ = 1.0 / ( liquidus_ - solidus_ ); + + temp_max_ = db.properties.vaporization_temperature; + rho_cp_ = rho * cp; rho_Lf_by_dT_ = rho * Lf / ( liquidus_ - solidus_ ); - k_by_dx2_ = ( db.properties.thermal_conductivity ) / ( dx * dx ); + inv_dx2_ = 1.0 / ( dx * dx ); + + k_solid_0_ = db.properties.thermal_conductivity_solid_0; + k_solid_T_ = db.properties.thermal_conductivity_solid_T; + k_liquid_0_ = db.properties.thermal_conductivity_liquid_0; + k_liquid_T_ = db.properties.thermal_conductivity_liquid_T; // initialize beam position for ( std::size_t d = 0; d < 3; ++d ) @@ -136,15 +151,52 @@ class Solver KOKKOS_INLINE_FUNCTION void operator()( HostTag tag, const int i, const int j, const int k ) const { - double x = T0_( i, j, k, 0 ); - - double dt_by_rho_cp = ( x >= solidus_ && x <= liquidus_ ) - ? dt_ / ( rho_cp_ + rho_Lf_by_dT_ ) - : dt_ / ( rho_cp_ ); + // First nearest neighbor stencil for cell at i,j,k: temperature and + // thermal conductivity + const double temp_local = T0_( i, j, k, 0 ); + const double temp_px = T0_( i + 1, j, k, 0 ); + const double temp_nx = T0_( i - 1, j, k, 0 ); + const double temp_py = T0_( i, j + 1, k, 0 ); + const double temp_ny = T0_( i, j - 1, k, 0 ); + const double temp_pz = T0_( i, j, k + 1, 0 ); + const double temp_nz = T0_( i, j, k - 1, 0 ); + + // Liquid fraction + const double liquid_fraction = + Kokkos::fmax( 1.0, Kokkos::fmin( 0.0, ( temp_local - solidus_ ) * + inv_freezing_range_ ) ); + const double kappa_local = + kappa_of_temperature( temp_local, liquid_fraction ); + const double kappa_px = + kappa_of_temperature( temp_px, liquid_fraction ); + const double kappa_nx = + kappa_of_temperature( temp_nx, liquid_fraction ); + const double kappa_py = + kappa_of_temperature( temp_py, liquid_fraction ); + const double kappa_ny = + kappa_of_temperature( temp_ny, liquid_fraction ); + const double kappa_pz = + kappa_of_temperature( temp_pz, liquid_fraction ); + const double kappa_nz = + kappa_of_temperature( temp_nz, liquid_fraction ); - double rhs = laplacian( i, j, k ) + source( tag, i, j, k ); - - T_( i, j, k, 0 ) = x + rhs * dt_by_rho_cp; + double dt_by_rho_cp = + ( liquid_fraction >= 0.0 && liquid_fraction <= 1.0 ) + ? dt_ / ( rho_cp_ + rho_Lf_by_dT_ ) + : dt_ / ( rho_cp_ ); + + const double laplacian_x = laplacian_k( + temp_local, temp_px, temp_nx, kappa_local, kappa_px, kappa_nx ); + const double laplacian_y = laplacian_k( + temp_local, temp_py, temp_ny, kappa_local, kappa_py, kappa_ny ); + const double laplacian_z = laplacian_k( + temp_local, temp_pz, temp_nz, kappa_local, kappa_pz, kappa_nz ); + + const double rhs = + ( laplacian_x + laplacian_y + laplacian_z ) * inv_dx2_ + + source( tag, i, j, k ); + + T_( i, j, k, 0 ) = temp_local + rhs * dt_by_rho_cp; } // Device tagged version of the temperature solver @@ -152,26 +204,85 @@ class Solver void operator()( DeviceTag tag, const int i, const int j, const int k ) const { - double x = T0_( i, j, k, 0 ); + // First nearest neighbor stencil for cell at i,j,k: temperature and + // thermal conductivity + const double temp_local = T0_( i, j, k, 0 ); + const double temp_px = T0_( i + 1, j, k, 0 ); + const double temp_nx = T0_( i - 1, j, k, 0 ); + const double temp_py = T0_( i, j + 1, k, 0 ); + const double temp_ny = T0_( i, j - 1, k, 0 ); + const double temp_pz = T0_( i, j, k + 1, 0 ); + const double temp_nz = T0_( i, j, k - 1, 0 ); + + // Liquid fraction + const double liquid_fraction = + Kokkos::fmax( 1.0, Kokkos::fmin( 0.0, ( temp_local - solidus_ ) * + inv_freezing_range_ ) ); + const double kappa_local = + kappa_of_temperature( temp_local, liquid_fraction ); + const double kappa_px = + kappa_of_temperature( temp_px, liquid_fraction ); + const double kappa_nx = + kappa_of_temperature( temp_nx, liquid_fraction ); + const double kappa_py = + kappa_of_temperature( temp_py, liquid_fraction ); + const double kappa_ny = + kappa_of_temperature( temp_ny, liquid_fraction ); + const double kappa_pz = + kappa_of_temperature( temp_pz, liquid_fraction ); + const double kappa_nz = + kappa_of_temperature( temp_nz, liquid_fraction ); double dt_by_rho_cp = - dt_ / ( rho_cp_ + - ( x >= solidus_ ) * ( x <= liquidus_ ) * rho_Lf_by_dT_ ); + dt_ / ( rho_cp_ + ( liquid_fraction >= 0.0 ) * + ( liquid_fraction <= 1.0 ) * rho_Lf_by_dT_ ); - double rhs = laplacian( i, j, k ) + source( tag, i, j, k ); + const double laplacian_x = laplacian_k( + temp_local, temp_px, temp_nx, kappa_local, kappa_px, kappa_nx ); + const double laplacian_y = laplacian_k( + temp_local, temp_py, temp_ny, kappa_local, kappa_py, kappa_ny ); + const double laplacian_z = laplacian_k( + temp_local, temp_pz, temp_nz, kappa_local, kappa_pz, kappa_nz ); - T_( i, j, k, 0 ) = x + rhs * dt_by_rho_cp; + const double rhs = + ( laplacian_x + laplacian_y + laplacian_z ) * inv_dx2_ + + source( tag, i, j, k ); + + T_( i, j, k, 0 ) = temp_local + rhs * dt_by_rho_cp; + } + + // Get temperature-dependent thermal conductivity, using liquid and solid + // values. Conductivity is capped at the value where temp_local = temp_max_ + KOKKOS_INLINE_FUNCTION + auto kappa_of_temperature( const double temp_local, + const double liquid_fraction ) const + { + return ( 1.0 - liquid_fraction ) * + ( k_solid_0_ + k_solid_T_ * temp_local ) + + liquid_fraction * + ( k_liquid_0_ + + k_liquid_T_ * Kokkos::fmin( temp_local, temp_max_ ) ); + } + + // Get harmonic average of two values + KOKKOS_INLINE_FUNCTION + auto harmonic_mean( const double x1, const double x2 ) const + { + return 2.0 * x1 * x2 / ( x1 + x2 ); } - // First-order centered space laplacian stencil + // First-order centered space laplacian stencil for one direction - + // temperature-dependent thermal conductivity KOKKOS_INLINE_FUNCTION - auto laplacian( const int i, const int j, const int k ) const + auto laplacian_k( const double temp_local, const double temp_positive, + const double temp_negative, const double kappa_local, + const double kappa_positive, + const double kappa_negative ) const { - return ( T0_( i - 1, j, k, 0 ) + T0_( i + 1, j, k, 0 ) + - T0_( i, j - 1, k, 0 ) + T0_( i, j + 1, k, 0 ) + - T0_( i, j, k - 1, 0 ) + T0_( i, j, k + 1, 0 ) - - 6.0 * T0_( i, j, k, 0 ) ) * - k_by_dx2_; + return ( harmonic_mean( kappa_local, kappa_positive ) * + ( temp_positive - temp_local ) - + harmonic_mean( kappa_local, kappa_negative ) * + ( temp_local - temp_negative ) ); } // Normalized weight for the gaussian source term: x in exp(-x)