Skip to content

Temperature-dependent thermal conductivity #26

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
MattRolchigo wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Temperature-dependent thermal conductivity #26

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9b764c3
Option for temperature-dependent thermal conductivity
MattRolchigo Dec 19, 2024
b7f3dc6
Separate conductivity for liquid and solid phases
MattRolchigo Dec 19, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 49 additions & 10 deletions src/Finch_Inputs.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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" ) &&
Copy link
Copy Markdown
Contributor

@streeve streeve Jan 2, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are these inputs "Mist compatible"?

Copy link
Copy Markdown
Collaborator Author

@MattRolchigo MattRolchigo Jan 7, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to Steve, Mist supports having properties represented as polynomials in the Y = a + b*X format so I think we're good there

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"];
Expand Down
155 changes: 133 additions & 22 deletions src/Finch_Solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_;
Expand All @@ -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 )
Expand Down Expand Up @@ -136,42 +151,138 @@ 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
KOKKOS_INLINE_FUNCTION
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)
Expand Down