diff --git a/examples/thermomechanics/inputs/thermal_deformation_heat_transfer.json b/examples/thermomechanics/inputs/thermal_deformation_heat_transfer.json index c3b68c5d..85abe80c 100644 --- a/examples/thermomechanics/inputs/thermal_deformation_heat_transfer.json +++ b/examples/thermomechanics/inputs/thermal_deformation_heat_transfer.json @@ -10,7 +10,7 @@ "horizon" : {"value": 0.00615, "unit": "m"}, "final_time" : {"value": 8, "unit": "s"}, "timestep" : {"value": 4.2e-07, "unit": "s"}, - "thermal_subcycle_steps" : {"value": 5e+4}, + "thermal_subcycle_steps" : {"value": 4.5e+4}, "output_frequency" : {"value": 10}, "output_reference" : {"value": true} } diff --git a/src/CabanaPD_Input.hpp b/src/CabanaPD_Input.hpp index 3105f2ee..d531b754 100644 --- a/src/CabanaPD_Input.hpp +++ b/src/CabanaPD_Input.hpp @@ -45,27 +45,18 @@ class Inputs int m = std::floor( delta / dx ); inputs["m"]["value"] = m; - // Set timestep safety factor if it not set by user + // Set timestep safety factor if not set by user. if ( !inputs.contains( "timestep_safety_factor" ) ) inputs["timestep_safety_factor"]["value"] = 0.85; - // Set approximate bulk modulus if not set by user - // Used in timestep estimation + // Check one of bulk or elastic modulus set by user. if ( !inputs.contains( "bulk_modulus" ) ) { if ( !inputs.contains( "elastic_modulus" ) ) throw std::runtime_error( "Must input either bulk_modulus or " "elastic_modulus." ); - double E = inputs["elastic_modulus"]["value"]; - double nu = 0.25; - double K = E / ( 3 * ( 1 - 2 * nu ) ); - inputs["bulk_modulus"]["value"] = K; } - // Check critical time step - // This must be done after the values above are calculated - computeCriticalTimeStep(); - int num_steps = tf / dt; inputs["num_steps"]["value"] = num_steps; @@ -177,72 +168,103 @@ class Inputs log_err( std::cout, "Units for dx do not match system units." ); } - void computeCriticalTimeStep() + template + void computeCriticalTimeStep( [[maybe_unused]] const ForceModel model ) { // Reference: Silling & Askari, Computers & Structures 83(17–18) (2005): // 1526-1535. - // Compute particle volume + // Compute particle volume. double dx = inputs["dx"]["value"][0]; double dy = inputs["dx"]["value"][1]; double dz = inputs["dx"]["value"][2]; double v_p = dx * dy * dz; - // Initialize denominator's summation + // Initialize denominator's summation. double sum = 0; + double sum_ht = 0; + + // Estimate the bulk modulus if needed. + double K; + if ( inputs.contains( "bulk_modulus" ) ) + { + K = inputs["bulk_modulus"]["value"]; + } + else + { + double E = inputs["elastic_modulus"]["value"]; + // This is only exact for bond-based (PMB). + double nu = 0.25; + K = E / ( 3 * ( 1 - 2 * nu ) ); + } - // Run over the neighborhood of a point in the bulk of a body + // Run over the neighborhood of a point in the bulk of a body (at the + // origin). int m = inputs["m"]["value"]; - double rho = inputs["density"]["value"]; - double K = inputs["bulk_modulus"]["value"]; double delta = inputs["horizon"]["value"]; // FIXME: this is copied from the forces double c = 18.0 * K / ( pi * delta * delta * delta * delta ); for ( int i = -( m + 1 ); i < m + 2; i++ ) { - // x-component of bond + // x-component of bond. double xi_1 = i * dx; for ( int j = -( m + 1 ); j < m + 2; j++ ) { - // y-component of bond + // y-component of bond. double xi_2 = j * dy; for ( int k = -( m + 1 ); k < m + 2; k++ ) { - // z-component of bond + // z-component of bond. double xi_3 = k * dz; - // Bond length squared + // Bond length squared. double r2 = xi_1 * xi_1 + xi_2 * xi_2 + xi_3 * xi_3; - // Check if bond is no longer than delta + // Check if bond is no longer than delta. if ( r2 < delta * delta + 1e-10 ) { - // Check is bond is not 0 + // Check if bond is not 0 if ( r2 > 0 ) { - // Compute denominator - sum += v_p * c / std::sqrt( r2 ); + double xi = std::sqrt( r2 ); + + // Compute denominator for mechanics. + sum += v_p * c / xi; + + // Compute denominator for heat transfer. + if constexpr ( is_heat_transfer< + typename ForceModel:: + thermal_type>::value ) + { + double coeff = + model.microconductivity_function( xi ); + sum_ht += v_p * coeff / r2; + } } } } } } - double safety_factor = inputs["timestep_safety_factor"]["value"]; - double dt_crit = safety_factor * std::sqrt( 2 * rho / sum ); - double dt = inputs["timestep"]["value"]; - if ( dt > dt_crit ) + double rho = inputs["density"]["value"]; + double dt_crit = std::sqrt( 2.0 * rho / sum ); + compareCriticalTimeStep( "mechanics", dt, dt_crit ); + + // Heat transfer timestep. + if constexpr ( is_heat_transfer< + typename ForceModel::thermal_type>::value ) { - log( std::cout, "WARNING: timestep (", dt, - ") is larger than estimated stable timestep (", dt_crit, - "), using safety factor of ", safety_factor, ".\n" ); + double dt_ht = inputs["thermal_subcycle_steps"]["value"]; + dt_ht *= dt; + + double cp = inputs["specific_heat_capacity"]["value"]; + double dt_ht_crit = rho * cp / sum_ht; + compareCriticalTimeStep( "heat_transfer", dt_ht, dt_ht_crit ); } - // Store in inputs - inputs["critical_timestep"]["value"] = dt_crit; } // Parse JSON file. @@ -268,6 +290,22 @@ class Inputs bool contains( std::string label ) { return inputs.contains( label ); } protected: + void compareCriticalTimeStep( std::string name, double dt, double dt_crit ) + { + double safety_factor = inputs["timestep_safety_factor"]["value"]; + double dt_crit_safety = safety_factor * dt_crit; + + if ( dt > dt_crit_safety ) + { + log( std::cout, "WARNING: timestep (", dt, + ") is larger than estimated stable timestep for ", name, " (", + dt_crit_safety, "), using safety factor of ", safety_factor, + ".\n" ); + } + // Store in inputs + inputs["critical_timestep_" + name]["value"] = dt_crit_safety; + } + nlohmann::json inputs; }; diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp index af642f53..50f7264f 100644 --- a/src/CabanaPD_Solver.hpp +++ b/src/CabanaPD_Solver.hpp @@ -142,6 +142,8 @@ class SolverElastic void setup( force_model_type force_model ) { + inputs.computeCriticalTimeStep( force_model ); + num_steps = inputs["num_steps"]; output_frequency = inputs["output_frequency"]; output_reference = inputs["output_reference"];