LCOV - code coverage report
Current view: top level - src/materials - ADGaussianHeatSourceBase.C (source / functions) Hit Total Coverage
Test: idaholab/malamute: 0e4c8a Lines: 100 118 84.7 %
Date: 2025-08-02 07:01:39 Functions: 7 8 87.5 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /****************************************************************************/
       2             : /*                        DO NOT MODIFY THIS HEADER                         */
       3             : /*                                                                          */
       4             : /* MALAMUTE: MOOSE Application Library for Advanced Manufacturing UTilitiEs */
       5             : /*                                                                          */
       6             : /*           Copyright 2021 - 2024, Battelle Energy Alliance, LLC           */
       7             : /*                           ALL RIGHTS RESERVED                            */
       8             : /****************************************************************************/
       9             : 
      10             : #include "ADGaussianHeatSourceBase.h"
      11             : 
      12             : InputParameters
      13          36 : ADGaussianHeatSourceBase::validParams()
      14             : {
      15          36 :   InputParameters params = Material::validParams();
      16          72 :   params.addRequiredParam<Real>("power", "laser power (1e-3 W)");
      17          72 :   params.addParam<Real>("efficiency", 1, "process efficiency");
      18          72 :   params.addParam<bool>("use_input_r",
      19          72 :                         true,
      20             :                         "option to use user input effective radii or from experimentally fitted "
      21             :                         "formulations. Default is to use user input data.");
      22          72 :   params.addParam<std::vector<Real>>("r",
      23             :                                      {},
      24             :                                      "effective radii (mm) along three directions. If only one "
      25             :                                      "parameter is provided, then we assume "
      26             :                                      "the effective radius to be equal along three directions.");
      27             : 
      28          72 :   params.addParam<Real>(
      29             :       "feed_rate",
      30          72 :       0.0,
      31             :       "powder material feed rate (g/ms). This value is used only when use_input_r = false.");
      32             : 
      33          72 :   params.addParam<Real>(
      34             :       "factor",
      35          72 :       1.0,
      36             :       "scaling factor that is multiplied to the heat source to adjust the intensity");
      37             : 
      38          72 :   params.addParam<Real>("std_factor", 1.0, "factor to the std to sample r");
      39             : 
      40          72 :   params.addParam<MooseEnum>(
      41         108 :       "heat_source_type", MooseEnum("point line mixed", "point"), "Type of the heat source");
      42             : 
      43          72 :   params.addParam<Real>("threshold_length",
      44          72 :                         1.0,
      45             :                         "Threshold size (mm) when we change the way of computing heat source");
      46             : 
      47         108 :   params.addParam<Real>(
      48             :       "number_time_integration",
      49          72 :       10,
      50             :       "Number of points to do time integration for averaged heat source calculation");
      51             : 
      52          72 :   params.declareControllable("power efficiency r factor");
      53             : 
      54          36 :   params.addClassDescription("Gaussian volumetric source heat source base.");
      55             : 
      56          36 :   return params;
      57           0 : }
      58             : 
      59          27 : ADGaussianHeatSourceBase::ADGaussianHeatSourceBase(const InputParameters & parameters)
      60             :   : Material(parameters),
      61          27 :     _use_input_r(getParam<bool>("use_input_r")),
      62          54 :     _P(getParam<Real>("power")),
      63          54 :     _eta(getParam<Real>("efficiency")),
      64          54 :     _feed_rate(getParam<Real>("feed_rate")),
      65          54 :     _r(getParam<std::vector<Real>>("r")),
      66          54 :     _f(getParam<Real>("factor")),
      67          54 :     _std_factor(getParam<Real>("std_factor")),
      68          54 :     _heat_source_type(getParam<MooseEnum>("heat_source_type").getEnum<HeatSourceType>()),
      69          54 :     _threshold_length(getParam<Real>("threshold_length")),
      70          54 :     _number_time_integration(getParam<Real>("number_time_integration")),
      71          54 :     _volumetric_heat(declareADProperty<Real>("volumetric_heat"))
      72             : {
      73          27 :   if (_use_input_r)
      74             :   {
      75          21 :     if (_r.size() != 1 && _r.size() != 3)
      76           0 :       paramError("r", "The effective radii should have 1 or 3 components.");
      77             :     // make sure we have 3 equal components if only one parameter is provided
      78          21 :     if (_r.size() == 1)
      79             :     {
      80          15 :       _r.push_back(_r[0]);
      81          15 :       _r.push_back(_r[0]);
      82             :     }
      83             :   }
      84             :   else
      85             :   {
      86           6 :     _r.resize(3);
      87             :     // Since the LINE and MIXED types averages over _t-_dt ->_t,
      88             :     // _r info for the previous steps & sub_steps are needed
      89             :     // however, this info is not available, and sampling again would change the
      90             :     // history
      91           6 :     if (_heat_source_type != HeatSourceType::POINT)
      92           0 :       paramError("heat_source_type",
      93             :                  "We can only use the POINT heat source type when _use_input_r = false.");
      94             :   }
      95             : 
      96             :   // set to a small number to start
      97          27 :   _r_time_prev = -1.0e8;
      98          27 : }
      99             : 
     100             : void
     101      282800 : ADGaussianHeatSourceBase::computeQpProperties()
     102             : {
     103      282800 :   const Real & x = _q_point[_qp](0);
     104             :   const Real & y = _q_point[_qp](1);
     105             :   const Real & z = _q_point[_qp](2);
     106             : 
     107      282800 :   switch (_heat_source_type)
     108             :   {
     109       80800 :     case HeatSourceType::POINT:
     110       80800 :       _volumetric_heat[_qp] = computeHeatSourceAtTime(x, y, z, _t);
     111       80800 :       break;
     112           0 :     case HeatSourceType::LINE:
     113           0 :       _volumetric_heat[_qp] = computeAveragedHeatSource(x, y, z, _t - _dt, _t);
     114           0 :       break;
     115      202000 :     case HeatSourceType::MIXED:
     116      202000 :       _volumetric_heat[_qp] = computeMixedHeatSource(x, y, z, _t - _dt, _t);
     117      202000 :       break;
     118             :   }
     119      282800 : }
     120             : 
     121             : Real
     122      282800 : ADGaussianHeatSourceBase::computeHeatSourceAtTime(const Real x,
     123             :                                                   const Real y,
     124             :                                                   const Real z,
     125             :                                                   const Real time)
     126             : {
     127             :   // center of the heat source
     128             :   Real x_t, y_t, z_t;
     129      282800 :   computeHeatSourceCenterAtTime(x_t, y_t, z_t, time);
     130             : 
     131      282800 :   Real dist_x = -2.0 * std::pow(x - x_t, 2.0) / _r[0] / _r[0];
     132      282800 :   Real dist_y = -2.0 * std::pow(y - y_t, 2.0) / _r[1] / _r[1];
     133      282800 :   Real dist_z = -2.0 * std::pow(z - z_t, 2.0) / _r[2] / _r[2];
     134             : 
     135             :   // Gaussian point heat source
     136      282800 :   return 2.0 * _P * _eta * _f / libMesh::pi / _r[0] / _r[1] / _r[2] *
     137      282800 :          std::exp(dist_x + dist_y + dist_z);
     138             : }
     139             : 
     140             : Real
     141           0 : ADGaussianHeatSourceBase::computeAveragedHeatSource(
     142             :     const Real x, const Real y, const Real z, const Real time_begin, const Real time_end)
     143             : {
     144             :   mooseAssert(time_end > time_begin, "Begin time should be smaller than end time.");
     145             :   // Use 5 points as a starting point. Number of points will double if
     146             :   // integration is not accurate enough or exceeds _number_time_integration.
     147             :   unsigned int num_pts = 5;
     148             :   Real Q_integral = 0, Q_integral_old = 0;
     149             :   do
     150             :   {
     151             :     Q_integral_old = Q_integral;
     152           0 :     Real delta_t = (time_end - time_begin) / num_pts;
     153             :     Real t0 = time_begin;
     154           0 :     Real Q_begin = computeHeatSourceAtTime(x, y, z, time_begin);
     155             :     Real Q_end;
     156             :     Q_integral = 0;
     157           0 :     for (unsigned int i = 0; i < num_pts; ++i)
     158             :     {
     159           0 :       t0 += delta_t;
     160           0 :       Q_end = computeHeatSourceAtTime(x, y, z, t0);
     161             :       // compute integral of Q between t0 and t0 + delta_t
     162           0 :       Q_integral += (Q_begin + Q_end) * delta_t / 2.0;
     163             :       // update Q_begin
     164             :       Q_begin = Q_end;
     165             :     }
     166           0 :     num_pts *= 2;
     167             :     // limit to _number_time_integration pts to accelerate the simulation
     168           0 :     if (num_pts > _number_time_integration)
     169             :       break;
     170           0 :   } while (!MooseUtils::absoluteFuzzyEqual(Q_integral_old, Q_integral, 1e-4));
     171             : 
     172           0 :   return Q_integral / (time_end - time_begin);
     173             : }
     174             : 
     175             : Real
     176      202000 : ADGaussianHeatSourceBase::computeMixedHeatSource(
     177             :     const Real x, const Real y, const Real z, const Real time_begin, const Real time_end)
     178             : {
     179             :   mooseAssert(time_end > time_begin, "Begin time should be smaller than end time.");
     180             : 
     181             :   // position at time_begin
     182             :   Real x_t0, y_t0, z_t0;
     183      202000 :   computeHeatSourceCenterAtTime(x_t0, y_t0, z_t0, time_begin);
     184      202000 :   Point P_start = Point(x_t0, y_t0, z_t0);
     185             :   // position at time_end
     186             :   Real x_t, y_t, z_t;
     187      202000 :   computeHeatSourceCenterAtTime(x_t, y_t, z_t, time_end);
     188      202000 :   Point P_end = Point(x_t, y_t, z_t);
     189             : 
     190      202000 :   Real SE = (P_end - P_start).norm();
     191      202000 :   if (SE < _threshold_length)
     192      202000 :     return computeHeatSourceAtTime(x, y, z, time_end);
     193             :   else
     194           0 :     return computeAveragedHeatSource(x, y, z, time_begin, time_end);
     195             : }
     196             : 
     197             : void
     198       35350 : ADGaussianHeatSourceBase::computeProperties()
     199             : {
     200             :   // effective radii under current processing parameters
     201       35350 :   if (!_use_input_r)
     202             :   {
     203             :     // compute scanning speed at this time
     204       10100 :     computeHeatSourceMovingSpeedAtTime(_t);
     205             :     // compute the effective radii
     206       10100 :     computeEffectiveRadii(_t);
     207             :   }
     208             : 
     209       35350 :   Material::computeProperties();
     210       35350 : }
     211             : 
     212             : void
     213       10100 : ADGaussianHeatSourceBase::computeEffectiveRadii(const Real time)
     214             : {
     215             :   // we do not update _r if we do not proceed in time
     216       10100 :   if (time <= _r_time_prev)
     217       10060 :     return;
     218             :   else
     219          40 :     _r_time_prev = time;
     220             : 
     221             :   // get scaled laser power, scanning speed, and powder feed rate
     222          40 :   Real lp = _P / 1.0e-3 / 400.0;             // input in unit 1e-3 W
     223          40 :   Real ss = _scan_speed / 4.23333e-4 / 40.0; // input in mm/ms (1 ipm = 4.23333e-4 mm/ms)
     224          40 :   Real pf = _feed_rate / 0.031e-3 / 15.0;    // input in g/ms (1rpm = 0.031e-3 g/ms)
     225             : 
     226             :   // list of the variable values
     227          40 :   std::vector<Real> vals = {lp, ss, pf, lp * lp, ss * ss, pf * pf, lp * ss, lp * pf, ss * pf, 1.0};
     228             : 
     229             :   Real mean_rxy = 0.0, mean_rz = 0.0;
     230             :   Real std_r = 0.0;
     231          40 :   _r[2] = 0.0;
     232         440 :   for (unsigned int i = 0; i < vals.size(); ++i)
     233             :   {
     234             :     // compute the mean and standard deviation of the melt pool dimension (x and y dimensions)
     235         400 :     mean_rxy += _diameter_param[i].first * vals[i];
     236         400 :     std_r += _diameter_param[i].second * vals[i];
     237             :     // compute the material bead height (z dimension)
     238         400 :     mean_rz += _height_param[i] * vals[i];
     239             :   }
     240             : 
     241             :   // sample the r[0] and r[2] value
     242             :   // define a random number generator
     243          40 :   std::random_device rd{};
     244          40 :   std::mt19937 generator{rd()};
     245             :   std::normal_distribution<double> dist_xy(mean_rxy, std_r);
     246             :   std::normal_distribution<double> dist_z(mean_rz, std_r);
     247          40 :   _r[0] = dist_xy(generator);
     248          40 :   _r[2] = dist_z(generator);
     249             :   // make sure that the sampled values are within +-sigma range
     250          40 :   if (_r[0] > mean_rxy + _std_factor * std::abs(std_r))
     251          20 :     _r[0] = mean_rxy + _std_factor * std::abs(std_r);
     252          20 :   else if (_r[0] < mean_rxy - _std_factor * std::abs(std_r) || _r[0] <= 0.0)
     253          20 :     _r[0] = std::abs(mean_rxy - _std_factor * std::abs(std_r)); // make sure r is not negative
     254          40 :   _r[0] *= 0.5;                                                 // get the radius
     255          40 :   _r[1] = _r[0];
     256             : 
     257          40 :   if (_r[2] > mean_rz + _std_factor * std::abs(std_r))
     258          17 :     _r[2] = mean_rz + _std_factor * std::abs(std_r);
     259          23 :   else if (_r[2] < mean_rz - _std_factor * std::abs(std_r) || _r[0] <= 0.0)
     260          23 :     _r[2] = std::abs(mean_rz - _std_factor * std::abs(std_r)); // make sure r is not negative
     261             : }

Generated by: LCOV version 1.14