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 - 2025, 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 17 : _r[0] = mean_rxy + _std_factor * std::abs(std_r);
252 23 : else if (_r[0] < mean_rxy - _std_factor * std::abs(std_r) || _r[0] <= 0.0)
253 23 : _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 21 : _r[2] = mean_rz + _std_factor * std::abs(std_r);
259 19 : else if (_r[2] < mean_rz - _std_factor * std::abs(std_r) || _r[0] <= 0.0)
260 19 : _r[2] = std::abs(mean_rz - _std_factor * std::abs(std_r)); // make sure r is not negative
261 40 : }
|
