diff --git a/2D_wire/2D_wire.e b/2D_wire/2D_wire.e new file mode 100644 index 0000000..2935338 Binary files /dev/null and b/2D_wire/2D_wire.e differ diff --git a/Makefile b/Makefile index 74da3bc..17f5128 100644 --- a/Makefile +++ b/Makefile @@ -34,7 +34,7 @@ FLUID_PROPERTIES := no FSI := no FUNCTIONAL_EXPANSION_TOOLS := no GEOCHEMISTRY := no -HEAT_TRANSFER := no +HEAT_TRANSFER := yes LEVEL_SET := no MISC := no NAVIER_STOKES := no diff --git a/axisymmetric_coil/AForm_axi_coil.i b/axisymmetric_coil/AForm_axi_coil.i index a78d421..cc0ee6a 100644 --- a/axisymmetric_coil/AForm_axi_coil.i +++ b/axisymmetric_coil/AForm_axi_coil.i @@ -1,7 +1,7 @@ [Mesh] [coil] type = FileMeshGenerator - file = axi_coil.e + file = joule_heating.e [] [tmg] type = TiledMeshGenerator @@ -24,7 +24,7 @@ [Kernels] [laplacian_A] # d2A/dz2 = 0 type = ADDiffusion - variable = A + variable = A [] [body_force] # d2A/dz2 = J*mu_0 type = ADBodyForce @@ -105,7 +105,7 @@ type = ParsedFunction expression = 'J' symbol_names = 'J' - symbol_values = '50' + symbol_values = '1e8' [] [] diff --git a/axisymmetric_coil/AForm_axi_coil_joule_heating.i b/axisymmetric_coil/AForm_axi_coil_joule_heating.i new file mode 100644 index 0000000..8d9e141 --- /dev/null +++ b/axisymmetric_coil/AForm_axi_coil_joule_heating.i @@ -0,0 +1,218 @@ +[Mesh] + [coil] + type = FileMeshGenerator + file = axi_coil.e + [] + [tmg] + type = TiledMeshGenerator + input = coil + left_boundary = left + right_boundary = right + top_boundary = top + bottom_boundary = bottom + x_tiles = 1 + y_tiles = 1 #change to number of turns in coil + [] + coord_type = rz + rz_coord_axis = y + [] + + [Problem] + type = FEProblem + [] + + [Kernels] + [laplacian_A] # d2A/dz2 = 0 + type = ADDiffusion + variable = A + [] + [body_force] # d2A/dz2 = J*mu_0 + type = ADBodyForce + variable = A + block = '2' + value = '1.257e-6' # mu_0 + function = 'J' + [] + [joule_heating] + type = ADJouleHeating #calculates Q = rho*J^2 + block = 2 #only the wire undergoes Joule heating due to presence of current density + variable = T + resistivity = resistivity + current_density = J + [] + [heat_conduction_tdot] + type = ADHeatConductionTimeDerivative #implements density * specific_heat * dT/dt + variable = T + [] + [heat_conduction] + type = ADHeatConduction #implements Q = -k*del(T) + variable = T + [] + [] + + [BCs] + [A_edge] + type = DirichletBC + variable = A + boundary = 'top left bottom right' + value = 0 + [] + [temp_edge] + type = DirichletBC + variable = T + boundary = 'top left bottom right' + value = 293.0 + [] + [convective_dissipation] + type = ADConvectiveHeatFluxBC + variable = T + boundary = 'wire_edge' + T_infinity = 293.0 + heat_transfer_coefficient = 10 + [] + [radiative_dissipation] + type = ADFunctionRadiativeBC + variable = T + boundary = 'wire_edge' + emissivity_function = '0.03' #approximate emissivity of copper wire + [] + [] + + [AuxKernels] + [dAdx_aux] + type = VariableGradientComponent + variable = dAdx + component = x + gradient_variable = A + [] + [dAdy_aux] + type = VariableGradientComponent + variable = dAdy + component = y + gradient_variable = A + [] + [mag_B_aux] + type = ParsedAux + variable = mag_B + coupled_variables = 'dAdx dAdy' + expression = 'sqrt(dAdx^2 + dAdy^2)' + [] + [Bx_unit_aux] + type = ParsedAux + variable = Bx_unit + coupled_variables = 'dAdy mag_B' + expression = 'dAdy / mag_B' + [] + [By_unit_aux] + type = ParsedAux + variable = By_unit + coupled_variables = 'dAdx mag_B' + expression = '-dAdx / mag_B' + [] + [] + + [AuxVariables] + [dAdx] + family = MONOMIAL + [] + [dAdy] + family = MONOMIAL + [] + [mag_B] + family = MONOMIAL + [] + [Bx_unit] + family = MONOMIAL + [] + [By_unit] + family = MONOMIAL + [] + [] + + [Variables] + [A] + [] + [T] + initial_condition = 293.0 + [] + [] + + [Materials] + [k_air] + type = ADGenericConstantMaterial + prop_names = 'thermal_conductivity' + prop_values = '0.03' #air in W/(m K) + block = 1 + [] + [cp_air] + type = ADGenericConstantMaterial + prop_names = 'specific_heat' + prop_values = '1000' #air in J/(kg K) + block = 1 + [] + [rho_air] + type = ADGenericConstantMaterial + prop_names = 'density' + prop_values = '1.293' #air in kg/(m^3) + block = 1 + [] + [k_copper] + type = ADGenericConstantMaterial + prop_names = 'thermal_conductivity' + prop_values = '397.48' #copper in W/(m K) + block = 2 + [] + [cp_copper] + type = ADGenericConstantMaterial + prop_names = 'specific_heat' + prop_values = '385.0' #copper in J/(kg K) + block = 2 + [] + [rho_copper] + type = ADGenericConstantMaterial + prop_names = 'density' + prop_values = '8920.0' #copper in kg/(m^3) + block = 2 + [] + [resistivity_copper] + type = ADResistivity + temperature = T + reference_resistivity = 1.68e-8 + temperature_coefficient = 0.00386 + reference_temperature = 293.0 + block = 2 + [] + [] + + [Functions] + [J] + type = ParsedFunction + expression = 'J' + symbol_names = 'J' + symbol_values = '1e8' + [] + [] + + [Preconditioning] + [SMP] + type = SMP + full = true + [] + [] + + [Executioner] + type = Transient + steady_state_detection = true + steady_state_tolerance = 1e-5 + scheme = bdf2 + solve_type = NEWTON + petsc_options_iname = '-pc_type' + petsc_options_value = 'hypre' + dt = 100 + end_time = 2000 + automatic_scaling = true + [] + + [Outputs] + exodus = true + [] \ No newline at end of file diff --git a/axisymmetric_coil/axi_coil.e b/axisymmetric_coil/axi_coil.e index bf7358e..2935338 100644 Binary files a/axisymmetric_coil/axi_coil.e and b/axisymmetric_coil/axi_coil.e differ diff --git a/include/kernels/ADJouleHeating.h b/include/kernels/ADJouleHeating.h new file mode 100644 index 0000000..b757a12 --- /dev/null +++ b/include/kernels/ADJouleHeating.h @@ -0,0 +1,34 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "ADKernelValue.h" + +class Function; + +/** + * This kernel calculates the heat source term corresponding to joule heating, + * Q = J * E = rho * J^2. + */ +class ADJouleHeating : public ADKernelValue +{ +public: + static InputParameters validParams(); + + ADJouleHeating(const InputParameters & parameters); + +protected: + virtual ADReal precomputeQpResidual() override; + +private: + const Function & _current_density; + + const ADMaterialProperty & _resistivity; +}; \ No newline at end of file diff --git a/include/materials/Resistivity.h b/include/materials/Resistivity.h new file mode 100644 index 0000000..dfea922 --- /dev/null +++ b/include/materials/Resistivity.h @@ -0,0 +1,42 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +#include "Material.h" + +/** + * Calculates resistivity as a function of temperature. + * It is assumed that resistivity varies linearly with temperature. + */ +template +class ResistivityTempl : public Material +{ +public: + static InputParameters validParams(); + + ResistivityTempl(const InputParameters & parameters); + +protected: + virtual void computeQpProperties(); + +private: + const Real _ref_resis; + const Real _temp_coeff; + const Real _ref_temp; + const bool _has_temp; + const GenericVariableValue & _T; + + const std::string _base_name; + GenericMaterialProperty & _resistivity; + MaterialProperty & _dresistivity_dT; +}; + +typedef ResistivityTempl Resistivity; +typedef ResistivityTempl ADResistivity; diff --git a/src/kernels/ADJouleHeating.C b/src/kernels/ADJouleHeating.C new file mode 100644 index 0000000..1b09b8b --- /dev/null +++ b/src/kernels/ADJouleHeating.C @@ -0,0 +1,41 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "ADJouleHeating.h" +#include "Function.h" + +registerMooseObject("orpheusApp", ADJouleHeating); + +InputParameters +ADJouleHeating::validParams() +{ + InputParameters params = ADKernelValue::validParams(); + params.addParam("current_density", "1", "A function that describes the current density"); + params.addParam( + "resistivity", + 1.77e-8, + "Material property providing resistivity of the material."); + params.addClassDescription("Calculates the heat source term corresponding to Joule " + "heating, with Jacobian contributions calculated using the automatic " + "differentiation system."); + return params; +} + +ADJouleHeating::ADJouleHeating(const InputParameters & parameters) + : ADKernelValue(parameters), + _current_density(getFunction("current_density")), + _resistivity(getADMaterialProperty("resistivity")) +{ +} + +ADReal +ADJouleHeating::precomputeQpResidual() +{ + return -_resistivity[_qp] * Utility::pow<2>(_current_density.value(_t, _q_point[_qp])); +} diff --git a/src/materials/Resistivity.C b/src/materials/Resistivity.C new file mode 100644 index 0000000..6cc805f --- /dev/null +++ b/src/materials/Resistivity.C @@ -0,0 +1,68 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#include "Resistivity.h" +#include "libmesh/quadrature.h" +#include "libmesh/utility.h" + +registerMooseObject("orpheusApp", Resistivity); +registerMooseObject("orpheusApp", ADResistivity); + +template +InputParameters +ResistivityTempl::validParams() +{ + InputParameters params = Material::validParams(); + params.addCoupledVar("temperature", 293.0, "Coupled variable for temperature."); + params.addParam("base_name", "Material property base name."); + params.addParam("reference_resistivity", + 1.68e-8, + "Electrical resistivity of the material at reference temperature " + "(default is copper resistivity in ohm-m)."); + params.addParam("temperature_coefficient", + 0.00386, + "Coefficient for calculating dependence of resistivity on temperature " + "(with copper as the default material)."); + params.addParam("reference_temperature", + 293.0, + "Reference temperature for Electrical resistivity in Kelvin."); + params.addClassDescription("Calculates resistivity as a function of " + "temperature, using copper for parameter defaults."); + return params; +} + +template +ResistivityTempl::ResistivityTempl(const InputParameters & parameters) + : Material(parameters), + _ref_resis(getParam("reference_resistivity")), + _temp_coeff(getParam("temperature_coefficient")), + _ref_temp(getParam("reference_temperature")), + _has_temp(isCoupled("temperature")), + _T(_has_temp ? coupledGenericValue("temperature") : genericZeroValue()), + _base_name(isParamValid("base_name") ? getParam("base_name") + "_" : ""), + _resistivity( + declareGenericProperty(_base_name + "resistivity")), + _dresistivity_dT(declareProperty(_base_name + "resistivity_dT")) +{ +} + +template +void +ResistivityTempl::computeQpProperties() +{ + const auto & temp_qp = _T[_qp]; + + const auto resistivity = _ref_resis * (1.0 + _temp_coeff * (temp_qp - _ref_temp)); + const Real dresistivity_dT = _ref_resis * _temp_coeff; + _resistivity[_qp] = resistivity; + _dresistivity_dT[_qp] = dresistivity_dT; +} + +template class ResistivityTempl; +template class ResistivityTempl; diff --git a/test/tests/kernels/joule_heating/gold/joule_heating_out.e b/test/tests/kernels/joule_heating/gold/joule_heating_out.e new file mode 100644 index 0000000..b0624a4 Binary files /dev/null and b/test/tests/kernels/joule_heating/gold/joule_heating_out.e differ diff --git a/test/tests/kernels/joule_heating/joule_heating.i b/test/tests/kernels/joule_heating/joule_heating.i new file mode 100644 index 0000000..2b4ee8f --- /dev/null +++ b/test/tests/kernels/joule_heating/joule_heating.i @@ -0,0 +1,135 @@ +[Mesh] + [coil] + type = FileMeshGenerator + file = joule_heating_mesh.e + [] + coord_type = rz + rz_coord_axis = y + [] + + [Problem] + type = FEProblem + [] + + [Kernels] + [joule_heating] + type = ADJouleHeating #calculates Q = rho*J^2 + block = 2 #only the wire undergoes Joule heating due to presence of current density + variable = T + resistivity = resistivity + current_density = J + [] + [heat_conduction_tdot] + type = ADHeatConductionTimeDerivative #implements density * specific_heat * dT/dt + variable = T + [] + [heat_conduction] + type = ADHeatConduction #implements Q = -k*del(T) + variable = T + [] + [] + + [BCs] + [temp_edge] + type = DirichletBC + variable = T + boundary = 'top left bottom right' + value = 293.0 + [] + [convective_dissipation] + type = ADConvectiveHeatFluxBC + variable = T + boundary = 'wire_edge' + T_infinity = 293.0 + heat_transfer_coefficient = 10 + [] + [radiative_dissipation] + type = ADFunctionRadiativeBC + variable = T + boundary = 'wire_edge' + emissivity_function = '0.03' #approximate emissivity of copper wire + [] + [] + + [Variables] + [T] + initial_condition = 293.0 + [] + [] + + [Materials] + [k_air] + type = ADGenericConstantMaterial + prop_names = 'thermal_conductivity' + prop_values = '0.03' #air in W/(m K) + block = 1 #air block + [] + [cp_air] + type = ADGenericConstantMaterial + prop_names = 'specific_heat' + prop_values = '1000' #air in J/(kg K) + block = 1 + [] + [rho_air] + type = ADGenericConstantMaterial + prop_names = 'density' + prop_values = '1.293' #air in kg/(m^3) + block = 1 + [] + [k_copper] + type = ADGenericConstantMaterial + prop_names = 'thermal_conductivity' + prop_values = '397.48' #copper in W/(m K) + block = 2 #copper block + [] + [cp_copper] + type = ADGenericConstantMaterial + prop_names = 'specific_heat' + prop_values = '385.0' #copper in J/(kg K) + block = 2 + [] + [rho_copper] + type = ADGenericConstantMaterial + prop_names = 'density' + prop_values = '8920.0' #copper in kg/(m^3) + block = 2 + [] + [resistivity_copper] + type = ADResistivity + temperature = T + block = 2 + [] + [] + + [Functions] + [J] + type = ParsedFunction + expression = 'J' + symbol_names = 'J' + symbol_values = '1e8' #current density + [] + [] + + [Preconditioning] + [SMP] + type = SMP + full = true + [] + [] + + [Executioner] + type = Transient + steady_state_detection = true + steady_state_tolerance = 1e-5 + scheme = bdf2 + solve_type = NEWTON + petsc_options_iname = '-pc_type' + petsc_options_value = 'hypre' + dt = 5 + end_time = 100 + automatic_scaling = true + [] + + [Outputs] + exodus = true + [] \ No newline at end of file diff --git a/test/tests/kernels/joule_heating/joule_heating_mesh.e b/test/tests/kernels/joule_heating/joule_heating_mesh.e new file mode 100644 index 0000000..2935338 Binary files /dev/null and b/test/tests/kernels/joule_heating/joule_heating_mesh.e differ diff --git a/test/tests/kernels/joule_heating/tests b/test/tests/kernels/joule_heating/tests new file mode 100644 index 0000000..e293c25 --- /dev/null +++ b/test/tests/kernels/joule_heating/tests @@ -0,0 +1,8 @@ +[Tests] + [test] + type = 'Exodiff' + input = 'joule_heating.i' + exodiff = 'joule_heating_out.e' + requirement = 'The system shall be able to solve a 2D axisymmetric Joule heating problem.' + [] +[] diff --git a/test/tests/materials/resistivity/gold/resistivity_test_out.e b/test/tests/materials/resistivity/gold/resistivity_test_out.e new file mode 100644 index 0000000..040e273 Binary files /dev/null and b/test/tests/materials/resistivity/gold/resistivity_test_out.e differ diff --git a/test/tests/materials/resistivity/resistivity_test.i b/test/tests/materials/resistivity/resistivity_test.i new file mode 100644 index 0000000..0e0d9b0 --- /dev/null +++ b/test/tests/materials/resistivity/resistivity_test.i @@ -0,0 +1,39 @@ +# Test for resistivity material property with copper as default values + +[Mesh] + type = GeneratedMesh + dim = 2 + nx = 2 + ny = 2 + xmax = 1 + ymax = 1 + [] + + [Problem] + solve = false + [] + + [Variables] + [T] + [] + [] + + [Materials] + [resistivity] + type = Resistivity + temperature = T + reference_resistivity = 1.68e-8 + temperature_coefficient = 0.00386 + reference_temperature = 293.0 + output_properties = 'resistivity dresistivity_dT' + outputs = exodus + [] + [] + + [Executioner] + type = Steady + [] + + [Outputs] + exodus = true + [] \ No newline at end of file diff --git a/test/tests/materials/resistivity/tests b/test/tests/materials/resistivity/tests new file mode 100644 index 0000000..3bfe8a0 --- /dev/null +++ b/test/tests/materials/resistivity/tests @@ -0,0 +1,8 @@ +[Tests] + [test] + type = 'Exodiff' + input = 'resistivity_test.i' + exodiff = 'resistivity_test_out.e' + [] +[] +