diff --git a/modules/heat_conduction/include/bcs/GaussianWeldEnergyFluxBC.h b/modules/heat_conduction/include/bcs/GaussianWeldEnergyFluxBC.h new file mode 100644 index 000000000000..0df9ce435fe2 --- /dev/null +++ b/modules/heat_conduction/include/bcs/GaussianWeldEnergyFluxBC.h @@ -0,0 +1,30 @@ +/****************************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */ +/* */ +/* Copyright 2021 - 2023, Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/****************************************************************************/ + +#pragma once + +#include "ADIntegratedBC.h" + +class GaussianWeldEnergyFluxBC : public ADIntegratedBC +{ +public: + static InputParameters validParams(); + + GaussianWeldEnergyFluxBC(const InputParameters & params); + +protected: + virtual ADReal computeQpResidual() override; + + const Real _reff; + const Real _F0; + const Real _R; + const Function & _x_beam_coord; + const Function & _y_beam_coord; + const Function & _z_beam_coord; +}; diff --git a/modules/heat_conduction/include/bcs/RadiationEnergyFluxBC.h b/modules/heat_conduction/include/bcs/RadiationEnergyFluxBC.h new file mode 100644 index 000000000000..d774504a4eff --- /dev/null +++ b/modules/heat_conduction/include/bcs/RadiationEnergyFluxBC.h @@ -0,0 +1,27 @@ +/****************************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */ +/* */ +/* Copyright 2021 - 2023, Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/****************************************************************************/ + +#pragma once + +#include "ADIntegratedBC.h" + +class RadiationEnergyFluxBC : public ADIntegratedBC +{ +public: + static InputParameters validParams(); + + RadiationEnergyFluxBC(const InputParameters & parameters); + +protected: + virtual ADReal computeQpResidual() override; + + const ADMaterialProperty & _sb_constant; + const ADMaterialProperty & _absorptivity; + const Real _ff_temp; +}; diff --git a/modules/heat_conduction/src/bcs/GaussianWeldEnergyFluxBC.C b/modules/heat_conduction/src/bcs/GaussianWeldEnergyFluxBC.C new file mode 100644 index 000000000000..72d8af19d969 --- /dev/null +++ b/modules/heat_conduction/src/bcs/GaussianWeldEnergyFluxBC.C @@ -0,0 +1,49 @@ +/****************************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */ +/* */ +/* Copyright 2021 - 2023, Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/****************************************************************************/ + +#include "GaussianWeldEnergyFluxBC.h" +#include "Function.h" + +registerMooseObject("MalamuteApp", GaussianWeldEnergyFluxBC); + +InputParameters +GaussianWeldEnergyFluxBC::validParams() +{ + InputParameters params = ADIntegratedBC::validParams(); + params.addRequiredParam("reff", + "The effective radius describing the radial distribution of the " + "beam energy. This should be non-dimensional."); + params.addRequiredParam("F0", "The average heat flux of the laser"); + params.addRequiredParam("R", "The beam radius"); + params.addParam("x_beam_coord", 0, "The x coordinate of the center of the beam"); + params.addParam("y_beam_coord", 0, "The y coordinate of the center of the beam"); + params.addParam("z_beam_coord", 0, "The z coordinate of the center of the beam"); + return params; +} + +GaussianWeldEnergyFluxBC::GaussianWeldEnergyFluxBC(const InputParameters & params) + : ADIntegratedBC(params), + _reff(getParam("reff")), + _F0(getParam("F0")), + _R(getParam("R")), + _x_beam_coord(getFunction("x_beam_coord")), + _y_beam_coord(getFunction("y_beam_coord")), + _z_beam_coord(getFunction("z_beam_coord")) +{ +} + +ADReal +GaussianWeldEnergyFluxBC::computeQpResidual() +{ + RealVectorValue beam_coords{_x_beam_coord.value(_t, _q_point[_qp]), + _y_beam_coord.value(_t, _q_point[_qp]), + _z_beam_coord.value(_t, _q_point[_qp])}; + auto r = (_ad_q_points[_qp] - beam_coords).norm(); + return -_test[_i][_qp] * 2. * _reff * _F0 * std::exp(-_reff * r * r / (_R * _R)); +} diff --git a/modules/heat_conduction/src/bcs/RadiationEnergyFluxBC.C b/modules/heat_conduction/src/bcs/RadiationEnergyFluxBC.C new file mode 100644 index 000000000000..983b3a26c90e --- /dev/null +++ b/modules/heat_conduction/src/bcs/RadiationEnergyFluxBC.C @@ -0,0 +1,42 @@ +/****************************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */ +/* */ +/* Copyright 2021 - 2023, Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/****************************************************************************/ + +#include "RadiationEnergyFluxBC.h" + +registerMooseObject("MalamuteApp", RadiationEnergyFluxBC); + +InputParameters +RadiationEnergyFluxBC::validParams() +{ + InputParameters params = ADIntegratedBC::validParams(); + params.addClassDescription("Computes heat flux due to radiation"); + params.addParam( + "sb_constant", "sb_constant", "The stefan-boltzmann constant"); + params.addParam("absorptivity", "abs", "The absorptivity of the material"); + params.addRequiredParam("ff_temp", "The far field temperature"); + return params; +} + +RadiationEnergyFluxBC::RadiationEnergyFluxBC(const InputParameters & parameters) + : ADIntegratedBC(parameters), + _sb_constant(getADMaterialProperty("sb_constant")), + _absorptivity(getADMaterialProperty("absorptivity")), + _ff_temp(getParam("ff_temp")) +{ +} + +ADReal +RadiationEnergyFluxBC::computeQpResidual() +{ + auto u2 = _u[_qp] * _u[_qp]; + auto u4 = u2 * u2; + auto ff2 = _ff_temp * _ff_temp; + auto ff4 = ff2 * ff2; + return _test[_i][_qp] * _absorptivity[_qp] * _sb_constant[_qp] * (u4 - ff4); +} diff --git a/modules/navier_stokes/include/bcs/DisplaceBoundaryBC.h b/modules/navier_stokes/include/bcs/DisplaceBoundaryBC.h new file mode 100644 index 000000000000..f7bd4daacc26 --- /dev/null +++ b/modules/navier_stokes/include/bcs/DisplaceBoundaryBC.h @@ -0,0 +1,26 @@ +/****************************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */ +/* */ +/* Copyright 2021 - 2023, Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/****************************************************************************/ + +#pragma once + +#include "ADNodalBC.h" + +class DisplaceBoundaryBC : public ADNodalBC +{ +public: + static InputParameters validParams(); + + DisplaceBoundaryBC(const InputParameters & parameters); + +protected: + virtual ADReal computeQpResidual() override; + + const ADReal & _velocity; + const Real & _u_old; +}; diff --git a/modules/navier_stokes/include/bcs/VaporRecoilPressureMomentumFluxBC.h b/modules/navier_stokes/include/bcs/VaporRecoilPressureMomentumFluxBC.h new file mode 100644 index 000000000000..8145592f7436 --- /dev/null +++ b/modules/navier_stokes/include/bcs/VaporRecoilPressureMomentumFluxBC.h @@ -0,0 +1,25 @@ +/****************************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */ +/* */ +/* Copyright 2021 - 2023, Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/****************************************************************************/ + +#pragma once + +#include "ADVectorIntegratedBC.h" + +class VaporRecoilPressureMomentumFluxBC : public ADVectorIntegratedBC +{ +public: + static InputParameters validParams(); + + VaporRecoilPressureMomentumFluxBC(const InputParameters & parameters); + +protected: + virtual ADReal computeQpResidual() override; + + const ADMaterialProperty & _rc_pressure; +}; diff --git a/modules/navier_stokes/src/bcs/DisplaceBoundaryBC.C b/modules/navier_stokes/src/bcs/DisplaceBoundaryBC.C new file mode 100644 index 000000000000..2d7504e65bb5 --- /dev/null +++ b/modules/navier_stokes/src/bcs/DisplaceBoundaryBC.C @@ -0,0 +1,34 @@ +/****************************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */ +/* */ +/* Copyright 2021 - 2023, Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/****************************************************************************/ + +#include "DisplaceBoundaryBC.h" + +registerMooseObject("MalamuteApp", DisplaceBoundaryBC); + +InputParameters +DisplaceBoundaryBC::validParams() +{ + InputParameters params = ADNodalBC::validParams(); + params.addClassDescription("For displacing a boundary"); + params.addRequiredCoupledVar("velocity", "The velocity at which to displace"); + return params; +} + +DisplaceBoundaryBC::DisplaceBoundaryBC(const InputParameters & parameters) + : ADNodalBC(parameters), + _velocity(adCoupledNodalValue("velocity")), + _u_old(_var.nodalValueOld()) +{ +} + +ADReal +DisplaceBoundaryBC::computeQpResidual() +{ + return _u - (_u_old + this->_dt * _velocity); +} diff --git a/modules/navier_stokes/src/bcs/VaporRecoilPressureMomentumFluxBC.C b/modules/navier_stokes/src/bcs/VaporRecoilPressureMomentumFluxBC.C new file mode 100644 index 000000000000..baa992e791f9 --- /dev/null +++ b/modules/navier_stokes/src/bcs/VaporRecoilPressureMomentumFluxBC.C @@ -0,0 +1,37 @@ +/****************************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */ +/* */ +/* Copyright 2021 - 2023, Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/****************************************************************************/ + +#include "VaporRecoilPressureMomentumFluxBC.h" + +registerMooseObject("MalamuteApp", VaporRecoilPressureMomentumFluxBC); + +InputParameters +VaporRecoilPressureMomentumFluxBC::validParams() +{ + InputParameters params = ADVectorIntegratedBC::validParams(); + params.addClassDescription("Vapor recoil pressure momentum flux"); + params.addParam("rc_pressure_name", "rc_pressure", "The recoil pressure"); + params.addCoupledVar("temperature", "The temperature on which the recoil pressure depends"); + return params; +} + +VaporRecoilPressureMomentumFluxBC::VaporRecoilPressureMomentumFluxBC( + const InputParameters & parameters) + : ADVectorIntegratedBC(parameters), _rc_pressure(getADMaterialProperty("rc_pressure_name")) +{ +} + +ADReal +VaporRecoilPressureMomentumFluxBC::computeQpResidual() +{ + return _test[_i][_qp] * + ADRealVectorValue( + std::abs(_normals[_qp](0)), std::abs(_normals[_qp](1)), std::abs(_normals[_qp](2))) * + _rc_pressure[_qp]; +} diff --git a/modules/navier_stokes/test/include/bcs/CrazyKCPlantFits.h b/modules/navier_stokes/test/include/bcs/CrazyKCPlantFits.h new file mode 100644 index 000000000000..70549f7cb27f --- /dev/null +++ b/modules/navier_stokes/test/include/bcs/CrazyKCPlantFits.h @@ -0,0 +1,52 @@ +/****************************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */ +/* */ +/* Copyright 2021 - 2023, Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/****************************************************************************/ + +#pragma once + +#include "ADMaterial.h" + +/** + * A material that couples a material property + */ +class CrazyKCPlantFits : public ADMaterial +{ +public: + static InputParameters validParams(); + + CrazyKCPlantFits(const InputParameters & parameters); + +protected: + virtual void computeQpProperties(); + + const Real _c_mu0; + const Real _c_mu1; + const Real _c_mu2; + const Real _c_mu3; + const Real _Tmax; + const Real _Tl; + const Real _T90; + const Real _beta; + const Real _c_k0; + const Real _c_k1; + const Real _c_cp0; + const Real _c_cp1; + const Real _c_rho0; + const ADVariableValue & _temperature; + const ADVariableGradient & _grad_temperature; + ADMaterialProperty & _mu; + ADMaterialProperty & _k; + ADMaterialProperty & _cp; + ADMaterialProperty & _rho; + ADMaterialProperty & _grad_k; + + const Real _length_units_per_meter; + const Real _temperature_units_per_kelvin; + const Real _mass_units_per_kilogram; + const Real _time_units_per_second; +}; diff --git a/modules/navier_stokes/test/include/bcs/CrazyKCPlantFitsBoundary.h b/modules/navier_stokes/test/include/bcs/CrazyKCPlantFitsBoundary.h new file mode 100644 index 000000000000..4ddd608a6248 --- /dev/null +++ b/modules/navier_stokes/test/include/bcs/CrazyKCPlantFitsBoundary.h @@ -0,0 +1,59 @@ +/****************************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */ +/* */ +/* Copyright 2021 - 2023, Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/****************************************************************************/ + +#pragma once + +#include "ADMaterial.h" + +/** + * A material that couples a material property + */ +class CrazyKCPlantFitsBoundary : public ADMaterial +{ +public: + static InputParameters validParams(); + + CrazyKCPlantFitsBoundary(const InputParameters & parameters); + +protected: + virtual void computeQpProperties(); + + const Real _ap0; + const Real _ap1; + const Real _ap2; + const Real _ap3; + const Real _bp0; + const Real _bp1; + const Real _bp2; + const Real _bp3; + const Real _Tb; + const Real _Tbound1; + const Real _Tbound2; + + const ADVariableValue & _temperature; + const ADVariableGradient & _grad_temperature; + + ADMaterialProperty & _rc_pressure; + + const Real _alpha; + const Real _sigma0; + const Real _T0; + ADMaterialProperty & _surface_tension; + ADMaterialProperty & _grad_surface_tension; + const MooseArray & _ad_normals; + const MooseArray & _ad_curvatures; + ADMaterialProperty & _surface_term_curvature; + ADMaterialProperty & _surface_term_gradient1; + ADMaterialProperty & _surface_term_gradient2; + + const Real _length_units_per_meter; + const Real _temperature_units_per_kelvin; + const Real _mass_units_per_kilogram; + const Real _time_units_per_second; +}; diff --git a/modules/navier_stokes/test/src/bcs/CrazyKCPlantFits.C b/modules/navier_stokes/test/src/bcs/CrazyKCPlantFits.C new file mode 100644 index 000000000000..02ca1e75f251 --- /dev/null +++ b/modules/navier_stokes/test/src/bcs/CrazyKCPlantFits.C @@ -0,0 +1,118 @@ +/****************************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */ +/* */ +/* Copyright 2021 - 2023, Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/****************************************************************************/ + +#include "CrazyKCPlantFits.h" + +registerMooseObject("MalamuteApp", CrazyKCPlantFits); + +InputParameters +CrazyKCPlantFits::validParams() +{ + InputParameters params = ADMaterial::validParams(); + params.addParam("c_mu0", 0.15616, "mu0 coefficient"); + params.addParam("c_mu1", -3.3696e-5, "mu1 coefficient"); + params.addParam("c_mu2", 1.0191e-8, "mu2 coefficient"); + params.addParam("c_mu3", -1.0413e-12, "mu3 coefficient"); + params.addParam("Tmax", 4000, "The maximum temperature"); + params.addParam("Tl", 1623, "The liquidus temperature"); + params.addParam( + "T90", 1528, "The T90 temperature (I don't know what this means physically)"); + params.addParam("beta", 1e11, "beta coefficient"); + params.addParam("c_k0", 10.7143, "k0 coefficient"); + params.addParam("c_k1", 14.2857e-3, "k0 coefficient"); + params.addParam("c_cp0", 425.75, "cp0 coefficient"); + params.addParam("c_cp1", 170.833e-3, "cp1 coefficient"); + params.addParam("c_rho0", 7.9e3, "The constant density"); + params.addCoupledVar("temperature", 1., "The temperature"); + params.addParam( + "mu_name", "mu", "The name of the viscosity material property"); + params.addParam("k_name", "k", "The name of the thermal conductivity"); + params.addParam("cp_name", "cp", "The name of the thermal conductivity"); + params.addParam("rho_name", "rho", "The name of the thermal conductivity"); + params.addParam("length_unit_exponent", + 0, + "The exponent of the length unit. If working in milimeters for example, " + "this number should be -3"); + params.addParam( + "temperature_unit_exponent", + 0, + "The exponent of the temperature unit. If working in kili-Kelvin for example, " + "this number should be 3"); + params.addParam("mass_unit_exponent", + 0, + "The exponent of the mass unit. If working in miligrams for example, " + "this number should be -9"); + params.addParam("time_unit_exponent", + 0, + "The exponent of the time unit. If working in micro-seconds for example, " + "this number should be -6"); + return params; +} + +CrazyKCPlantFits::CrazyKCPlantFits(const InputParameters & parameters) + : ADMaterial(parameters), + _c_mu0(getParam("c_mu0")), + _c_mu1(getParam("c_mu1")), + _c_mu2(getParam("c_mu2")), + _c_mu3(getParam("c_mu3")), + _Tmax(getParam("Tmax")), + _Tl(getParam("Tl")), + _T90(getParam("T90")), + _beta(getParam("beta")), + _c_k0(getParam("c_k0")), + _c_k1(getParam("c_k1")), + _c_cp0(getParam("c_cp0")), + _c_cp1(getParam("c_cp1")), + _c_rho0(getParam("c_rho0")), + _temperature(adCoupledValue("temperature")), + _grad_temperature(adCoupledGradient("temperature")), + _mu(declareADProperty(getParam("mu_name"))), + _k(declareADProperty(getParam("k_name"))), + _cp(declareADProperty(getParam("cp_name"))), + _rho(declareADProperty(getParam("rho_name"))), + _grad_k(declareADProperty("grad_" + getParam("k_name"))), + _length_units_per_meter(1. / std::pow(10, getParam("length_unit_exponent"))), + _temperature_units_per_kelvin(1. / std::pow(10, getParam("temperature_unit_exponent"))), + _mass_units_per_kilogram(1. / std::pow(10, getParam("mass_unit_exponent"))), + _time_units_per_second(1. / std::pow(10, getParam("time_unit_exponent"))) +{ +} + +void +CrazyKCPlantFits::computeQpProperties() +{ + if (_temperature[_qp] < _Tl * _temperature_units_per_kelvin) + _mu[_qp] = _mass_units_per_kilogram / (_length_units_per_meter * _time_units_per_second) * + (_c_mu0 + _c_mu1 * _Tl + _c_mu2 * _Tl * _Tl + _c_mu3 * _Tl * _Tl * _Tl) * + (_beta + (1 - _beta) * (_temperature[_qp] - _T90 * _temperature_units_per_kelvin) / + ((_Tl - _T90) * _temperature_units_per_kelvin)); + else + { + ADReal That; + That = _temperature[_qp] / _temperature_units_per_kelvin > _Tmax + ? _Tmax + : _temperature[_qp] / _temperature_units_per_kelvin; + _mu[_qp] = _mass_units_per_kilogram / (_length_units_per_meter * _time_units_per_second) * + (_c_mu0 + _c_mu1 * That + _c_mu2 * That * That + _c_mu3 * That * That * That); + } + _k[_qp] = (_c_k0 + _c_k1 / _temperature_units_per_kelvin * _temperature[_qp]) * + (_mass_units_per_kilogram * _length_units_per_meter / + (_temperature_units_per_kelvin * _time_units_per_second * _time_units_per_second * + _time_units_per_second)); + _grad_k[_qp] = _c_k1 * + (_mass_units_per_kilogram * _length_units_per_meter / + (_temperature_units_per_kelvin * _temperature_units_per_kelvin * + _time_units_per_second * _time_units_per_second * _time_units_per_second)) * + _grad_temperature[_qp]; + _cp[_qp] = (_c_cp0 + _c_cp1 / _temperature_units_per_kelvin * _temperature[_qp]) * + (_length_units_per_meter * _length_units_per_meter / + (_temperature_units_per_kelvin * _time_units_per_second * _time_units_per_second)); + _rho[_qp] = _c_rho0 * _mass_units_per_kilogram / + (_length_units_per_meter * _length_units_per_meter * _length_units_per_meter); +} diff --git a/modules/navier_stokes/test/src/bcs/CrazyKCPlantFitsBoundary.C b/modules/navier_stokes/test/src/bcs/CrazyKCPlantFitsBoundary.C new file mode 100644 index 000000000000..776bf36bec71 --- /dev/null +++ b/modules/navier_stokes/test/src/bcs/CrazyKCPlantFitsBoundary.C @@ -0,0 +1,124 @@ +/****************************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */ +/* */ +/* Copyright 2021 - 2023, Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/****************************************************************************/ + +#include "CrazyKCPlantFitsBoundary.h" +#include "Assembly.h" + +registerMooseObject("MalamuteApp", CrazyKCPlantFitsBoundary); + +InputParameters +CrazyKCPlantFitsBoundary::validParams() +{ + InputParameters params = ADMaterial::validParams(); + params.addParam("c_mu0", 0.15616, "mu0 coefficient"); + params.addParam("ap0", 0, ""); + params.addParam("ap1", 1.851502e1, ""); + params.addParam("ap2", -1.96945e-1, ""); + params.addParam("ap3", 1.594124e-3, ""); + params.addParam("bp0", 0, ""); + params.addParam("bp1", -5.809553e1, ""); + params.addParam("bp2", 4.610515e-1, ""); + params.addParam("bp3", 2.332819e-4, ""); + params.addParam("Tb", 3000, "The boiling temperature"); + params.addParam("Tbound1", 0, "The first temperature bound"); + params.addParam("Tbound2", 170, "The second temperature bound"); + params.addCoupledVar("temperature", 1., "The temperature"); + params.addParam("rc_pressure_name", "rc_pressure", "The recoil pressure"); + params.addParam( + "alpha", -4.3e-4, "The derivative of the surface tension with respect to temperature"); + params.addParam("sigma0", 1.943, "The surface tension at T0"); + params.addParam("T0", 1809, "The reference temperature for the surface tension"); + params.addParam( + "surface_tension_name", "surface_tension", "The surface tension"); + params.addParam( + "grad_surface_tension_name", "grad_surface_tension", "The gradient of the surface tension"); + params.addParam("length_unit_exponent", + 0, + "The exponent of the length unit. If working in milimeters for example, " + "this number should be -3"); + params.addParam( + "temperature_unit_exponent", + 0, + "The exponent of the temperature unit. If working in kili-Kelvin for example, " + "this number should be 3"); + params.addParam("mass_unit_exponent", + 0, + "The exponent of the mass unit. If working in miligrams for example, " + "this number should be -9"); + params.addParam("time_unit_exponent", + 0, + "The exponent of the time unit. If working in micro-seconds for example, " + "this number should be -6"); + return params; +} + +CrazyKCPlantFitsBoundary::CrazyKCPlantFitsBoundary(const InputParameters & parameters) + : ADMaterial(parameters), + _ap0(getParam("ap0")), + _ap1(getParam("ap1")), + _ap2(getParam("ap2")), + _ap3(getParam("ap3")), + _bp0(getParam("bp0")), + _bp1(getParam("bp1")), + _bp2(getParam("bp2")), + _bp3(getParam("bp3")), + _Tb(getParam("Tb")), + _Tbound1(getParam("Tbound1")), + _Tbound2(getParam("Tbound2")), + _temperature(adCoupledValue("temperature")), + _grad_temperature(adCoupledGradient("temperature")), + _rc_pressure(declareADProperty(getParam("rc_pressure_name"))), + _alpha(getParam("alpha")), + _sigma0(getParam("sigma0")), + _T0(getParam("T0")), + _surface_tension( + declareADProperty(getParam("surface_tension_name"))), + _grad_surface_tension(declareADProperty( + getParam("grad_surface_tension_name"))), + _ad_normals(_assembly.adNormals()), + _ad_curvatures(_assembly.adCurvatures()), + _surface_term_curvature(declareADProperty("surface_term_curvature")), + _surface_term_gradient1(declareADProperty("surface_term_gradient1")), + _surface_term_gradient2(declareADProperty("surface_term_gradient2")), + _length_units_per_meter(1. / std::pow(10, getParam("length_unit_exponent"))), + _temperature_units_per_kelvin(1. / std::pow(10, getParam("temperature_unit_exponent"))), + _mass_units_per_kilogram(1. / std::pow(10, getParam("mass_unit_exponent"))), + _time_units_per_second(1. / std::pow(10, getParam("time_unit_exponent"))) +{ +} + +void +CrazyKCPlantFitsBoundary::computeQpProperties() +{ + auto && theta = _temperature[_qp] / _temperature_units_per_kelvin - _Tb; + if (theta < _Tbound1) + _rc_pressure[_qp] = 0; + else if (theta < _Tbound2) + _rc_pressure[_qp] = + _mass_units_per_kilogram / + (_length_units_per_meter * _time_units_per_second * _time_units_per_second) * + (_ap0 + _ap1 * theta + _ap2 * theta * theta + _ap3 * theta * theta * theta); + else + _rc_pressure[_qp] = + _mass_units_per_kilogram / + (_length_units_per_meter * _time_units_per_second * _time_units_per_second) * + (_bp0 + _bp1 * theta + _bp2 * theta * theta + _bp3 * theta * theta * theta); + + _surface_tension[_qp] = + _sigma0 * _mass_units_per_kilogram / (_time_units_per_second * _time_units_per_second) + + _alpha * _mass_units_per_kilogram / (_time_units_per_second * _time_units_per_second) * + (_temperature[_qp] / _temperature_units_per_kelvin - _T0); + _grad_surface_tension[_qp] = _alpha * _mass_units_per_kilogram / + (_time_units_per_second * _time_units_per_second) / + _temperature_units_per_kelvin * _grad_temperature[_qp]; + _surface_term_curvature[_qp] = + -2. * _ad_curvatures[_qp] * _surface_tension[_qp] * _ad_normals[_qp]; + _surface_term_gradient1[_qp] = -_grad_surface_tension[_qp]; + _surface_term_gradient2[_qp] = _ad_normals[_qp] * (_ad_normals[_qp] * _grad_surface_tension[_qp]); +} diff --git a/modules/navier_stokes/test/tests/finite_element/ins/laser-welding/3d.i b/modules/navier_stokes/test/tests/finite_element/ins/laser-welding/3d.i new file mode 100644 index 000000000000..339146c01dec --- /dev/null +++ b/modules/navier_stokes/test/tests/finite_element/ins/laser-welding/3d.i @@ -0,0 +1,511 @@ +period=1.25e-3 +endtime=2.5e-3 +timestep=1.25e-5 +surfacetemp=300 + +[GlobalParams] + gravity = '0 0 0' + pspg = true + # supg = true + laplace = true + integrate_p_by_parts = true + convective_term = true + transient_term = true + temperature = T +[] + +[Mesh] + type = GeneratedMesh + dim = 3 + xmin = -.35e-3 + xmax = 0.35e-3 + ymin = -.35e-3 + ymax = .35e-3 + zmin = -.7e-3 + zmax = 0 + nx = 4 + ny = 4 + nz = 4 + displacements = 'disp_x disp_y disp_z' + uniform_refine = 2 +[] + +[Variables] + [./vel_x] + [./InitialCondition] + type = ConstantIC + value = 1e-15 + [../] + [../] + + [./vel_y] + [./InitialCondition] + type = ConstantIC + value = 1e-15 + [../] + [../] + + [./vel_z] + [./InitialCondition] + type = ConstantIC + value = 1e-15 + [../] + [../] + + [./T] + [../] + + [./p] + [../] + [./disp_x] + [../] + [./disp_y] + [../] + [./disp_z] + [../] +[] + +[ICs] + [./T] + type = FunctionIC + variable = T + function = '(${surfacetemp} - 300) / .7e-3 * z + ${surfacetemp}' + [../] +[] + +[Kernels] + [./disp_x] + type = Diffusion + variable = disp_x + [../] + [./disp_y] + type = Diffusion + variable = disp_y + [../] + [./disp_z] + type = Diffusion + variable = disp_z + [../] +[] + +[ADKernels] + [./mesh_x] + type = INSConvectedMesh + variable = vel_x + disp_x = disp_x + disp_y = disp_y + disp_z = disp_z + use_displaced_mesh = true + [../] + [./mesh_y] + type = INSConvectedMesh + variable = vel_y + disp_x = disp_x + disp_y = disp_y + disp_z = disp_z + use_displaced_mesh = true + [../] + [./mesh_z] + type = INSConvectedMesh + variable = vel_z + disp_x = disp_x + disp_y = disp_y + disp_z = disp_z + use_displaced_mesh = true + [../] + + # mass + [./mass] + type = INSADMass + variable = p + u = vel_x + v = vel_y + w = vel_z + p = p + use_displaced_mesh = true + [../] + + # x-momentum, time + [./x_momentum_time] + type = INSADMomentumTimeDerivative + variable = vel_x + use_displaced_mesh = true + [../] + + # x-momentum, space + [./x_momentum_space] + type = INSADMomentumBase + variable = vel_x + u = vel_x + v = vel_y + w = vel_z + p = p + component = 0 + use_displaced_mesh = true + [../] + + # y-momentum, time + [./y_momentum_time] + type = INSADMomentumTimeDerivative + variable = vel_y + use_displaced_mesh = true + [../] + + # y-momentum, space + [./y_momentum_space] + type = INSADMomentumBase + variable = vel_y + u = vel_x + v = vel_y + w = vel_z + p = p + component = 1 + use_displaced_mesh = true + [../] + + # z-momentum, time + [./z_momentum_time] + type = INSADMomentumTimeDerivative + variable = vel_z + use_displaced_mesh = true + [../] + + # z-momentum, space + [./z_momentum_space] + type = INSADMomentumBase + variable = vel_z + u = vel_x + v = vel_y + w = vel_z + p = p + component = 2 + use_displaced_mesh = true + [../] + + # temperature + [./temperature_time] + type = INSADTemperatureTimeDerivative + variable = T + use_displaced_mesh = true + [../] + + [./temperature_space] + type = INSADTemperature + variable = T + u = vel_x + v = vel_y + w = vel_z + p = p + use_displaced_mesh = true + [../] + [mesh_T] + type = INSTemperatureConvectedMesh + variable = T + disp_x = disp_x + disp_y = disp_y + disp_z = disp_z + use_displaced_mesh = true + [] +[] + +[BCs] + [./x_no_disp] + type = DirichletBC + variable = disp_x + boundary = 'back' + value = 0 + [../] + [./y_no_disp] + type = DirichletBC + variable = disp_y + boundary = 'back' + value = 0 + [../] + [./z_no_disp] + type = DirichletBC + variable = disp_z + boundary = 'back' + value = 0 + [../] + + [./x_no_slip] + type = DirichletBC + variable = vel_x + boundary = 'bottom right left top back' + value = 0.0 + [../] + + [./y_no_slip] + type = DirichletBC + variable = vel_y + boundary = 'bottom right left top back' + value = 0.0 + [../] + + [./z_no_slip] + type = DirichletBC + variable = vel_z + boundary = 'bottom right left top back' + value = 0.0 + [../] + + [./T_cold] + type = DirichletBC + variable = T + boundary = 'back' + value = 300 + [../] +[] + +[ADBCs] + [./radiation_flux] + type = RadiationEnergyFluxBC + variable = T + boundary = 'front' + ff_temp = 300 + sb_constant = 'sb_constant' + absorptivity = 'abs' + use_displaced_mesh = true + [../] + + [./weld_flux] + type = GaussianWeldEnergyFluxBC + variable = T + boundary = 'front' + reff = 0.6 + F0 = 2.546e9 + R = 1e-4 + x_beam_coord = '2e-4 * cos(t * 2 * pi / ${period})' + y_beam_coord = '2e-4 * sin(t * 2 * pi / ${period})' + z_beam_coord = 0 + use_displaced_mesh = true + [../] + + [./vapor_recoil_x] + type = VaporRecoilPressureMomentumFluxBC + variable = vel_x + boundary = 'front' + component = 0 + use_displaced_mesh = true + [../] + + [./vapor_recoil_y] + type = VaporRecoilPressureMomentumFluxBC + variable = vel_y + boundary = 'front' + component = 1 + use_displaced_mesh = true + [../] + + [./vapor_recoil_z] + type = VaporRecoilPressureMomentumFluxBC + variable = vel_z + boundary = 'front' + component = 2 + use_displaced_mesh = true + [../] + +[./displace_x_top] + type = DisplaceBoundaryBC + boundary = 'front' + variable = 'disp_x' + velocity = 'vel_x' + [../] + [./displace_y_top] + type = DisplaceBoundaryBC + boundary = 'front' + variable = 'disp_y' + velocity = 'vel_y' + [../] + [./displace_z_top] + type = DisplaceBoundaryBC + boundary = 'front' + variable = 'disp_z' + velocity = 'vel_z' + [../] +[] + +[ADMaterials] + [./kc_fits] + type = CrazyKCPlantFits + temperature = T + beta = 1e7 + [../] + [./boundary] + type = CrazyKCPlantFitsBoundary + boundary = 'front' + temperature = T + [../] +[] + +[Materials] + [./const] + type = GenericConstantMaterial + prop_names = 'abs sb_constant' + prop_values = '1 5.67e-8' + [../] +[] + +[Preconditioning] + [./SMP] + type = SMP + full = true + solve_type = 'NEWTON' + [../] +[] + +[Executioner] + type = Transient + end_time = ${endtime} + dtmin = 1e-8 + dtmax = ${timestep} + petsc_options = '-snes_converged_reason -ksp_converged_reason -options_left -ksp_monitor_singular_value' + petsc_options_iname = '-ksp_gmres_restart -pc_type' + petsc_options_value = '100 asm' + + line_search = 'none' + nl_max_its = 12 + l_max_its = 100 + [./TimeStepper] + type = IterationAdaptiveDT + optimal_iterations = 7 + dt = ${timestep} + linear_iteration_ratio = 1e6 + growth_factor = 1.5 + [../] +[] + +[Outputs] + print_linear_residuals = false + [./exodus] + type = Exodus + output_material_properties = true + show_material_properties = 'mu' + [../] + [./dofmap] + type = DOFMap + execute_on = 'initial' + [../] + checkpoint = true +[] + +[Debug] + show_var_residual_norms = true +[] + + +[Adaptivity] + marker = combo + max_h_level = 4 + + [./Indicators] + # [./error_x] + # type = GradientJumpIndicator + # variable = vel_x + # [../] + # [./error_y] + # type = GradientJumpIndicator + # variable = vel_y + # [../] + # [./error_z] + # type = GradientJumpIndicator + # variable = vel_z + # [../] + # [./error_p] + # type = GradientJumpIndicator + # variable = p + # [../] + [./error_T] + type = GradientJumpIndicator + variable = T + [../] + # [./error_dispx] + # type = GradientJumpIndicator + # variable = disp_x + # [../] + # [./error_dispy] + # type = GradientJumpIndicator + # variable = disp_y + # [../] + [./error_dispz] + type = GradientJumpIndicator + variable = disp_z + [../] + [../] + + [./Markers] + # [./errorfrac_x] + # type = ErrorFractionMarker + # refine = 0.7 + # coarsen = 0.3 + # indicator = error_x + # [../] + # [./errorfrac_y] + # type = ErrorFractionMarker + # refine = 0.7 + # coarsen = 0.3 + # indicator = error_y + # [../] + # [./errorfrac_z] + # type = ErrorFractionMarker + # refine = 0.7 + # coarsen = 0.3 + # indicator = error_z + # [../] + # [./errorfrac_p] + # type = ErrorFractionMarker + # refine = 0.7 + # coarsen = 0.3 + # indicator = error_p + # [../] + [./errorfrac_T] + type = ErrorFractionMarker + refine = 0.4 + coarsen = 0.2 + indicator = error_T + [../] + # [./errorfrac_dispx] + # type = ErrorFractionMarker + # refine = 0.7 + # coarsen = 0.3 + # indicator = error_dispx + # [../] + # [./errorfrac_dispy] + # type = ErrorFractionMarker + # refine = 0.7 + # coarsen = 0.3 + # indicator = error_dispy + # [../] + [./errorfrac_dispz] + type = ErrorFractionMarker + refine = 0.4 + coarsen = 0.2 + indicator = error_dispz + [../] + [./combo] + type = ComboMarker + markers = 'errorfrac_T errorfrac_dispz' + [../] + [../] +[] + +[Postprocessors] + [./num_dofs] + type = NumDOFs + system = 'NL' + [../] + [./nl] + type = NumNonlinearIterations + [../] + [./tot_nl] + type = CumulativeValuePostprocessor + postprocessor = 'nl' + [../] + [./linear] + type = NumLinearIterations + [../] + [./tot_linear] + type = CumulativeValuePostprocessor + postprocessor = 'linear' + [../] +[]