Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "TriInterWrapper1PhaseProblem.h"
11 : #include "AuxiliarySystem.h"
12 : #include "TriInterWrapperMesh.h"
13 : #include "SCM.h"
14 :
15 : registerMooseObject("SubChannelApp", TriInterWrapper1PhaseProblem);
16 :
17 : InputParameters
18 78 : TriInterWrapper1PhaseProblem::validParams()
19 : {
20 78 : InputParameters params = InterWrapper1PhaseProblem::validParams();
21 78 : params.addClassDescription(
22 : "Solver class for interwrapper of assemblies in a triangular-lattice arrangement");
23 78 : return params;
24 0 : }
25 :
26 39 : TriInterWrapper1PhaseProblem::TriInterWrapper1PhaseProblem(const InputParameters & params)
27 : : InterWrapper1PhaseProblem(params),
28 39 : _tri_sch_mesh(SCM::getConstMesh<TriInterWrapperMesh>(_subchannel_mesh))
29 : {
30 39 : }
31 :
32 : double
33 504000 : TriInterWrapper1PhaseProblem::computeFrictionFactor(Real Re)
34 : {
35 : Real a, b;
36 504000 : if (Re < 1)
37 : {
38 : return 64.0;
39 : }
40 504000 : else if (Re >= 1 and Re < 5000)
41 : {
42 : a = 64.0;
43 : b = -1.0;
44 : }
45 0 : else if (Re >= 5000 and Re < 30000)
46 : {
47 : a = 0.316;
48 : b = -0.25;
49 : }
50 : else
51 : {
52 : a = 0.184;
53 : b = -0.20;
54 : }
55 504000 : return a * std::pow(Re, b);
56 : }
57 :
58 : void
59 1494 : TriInterWrapper1PhaseProblem::computeDP(int iblock)
60 : {
61 1494 : unsigned int last_node = (iblock + 1) * _block_size;
62 1494 : unsigned int first_node = iblock * _block_size + 1;
63 :
64 13494 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
65 : {
66 12000 : auto z_grid = _subchannel_mesh.getZGrid();
67 12000 : auto k_grid = _subchannel_mesh.getKGrid();
68 12000 : auto dz = z_grid[iz] - z_grid[iz - 1];
69 516000 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
70 : {
71 504000 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
72 504000 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
73 504000 : auto rho_in = (*_rho_soln)(node_in);
74 504000 : auto rho_out = (*_rho_soln)(node_out);
75 504000 : auto mu_in = (*_mu_soln)(node_in);
76 504000 : auto S = (*_S_flow_soln)(node_in);
77 504000 : auto w_perim = (*_w_perim_soln)(node_in);
78 : // hydraulic diameter in the i direction
79 504000 : auto Dh_i = 4.0 * S / w_perim;
80 : auto time_term =
81 504000 : _TR * ((*_mdot_soln)(node_out)-_mdot_soln->old(node_out)) * dz / _dt -
82 504000 : dz * 2.0 * (*_mdot_soln)(node_out) * (rho_out - _rho_soln->old(node_out)) / rho_in / _dt;
83 :
84 : auto mass_term1 =
85 504000 : std::pow((*_mdot_soln)(node_out), 2.0) * (1.0 / S / rho_out - 1.0 / S / rho_in);
86 504000 : auto mass_term2 = -2.0 * (*_mdot_soln)(node_out) * (*_SumWij_soln)(node_out) / S / rho_in;
87 :
88 : auto crossflow_term = 0.0;
89 : auto turbulent_term = 0.0;
90 :
91 : unsigned int counter = 0;
92 1944000 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
93 : {
94 1440000 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
95 : unsigned int ii_ch = chans.first;
96 : unsigned int jj_ch = chans.second;
97 1440000 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
98 1440000 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
99 1440000 : auto * node_out_i = _subchannel_mesh.getChannelNode(ii_ch, iz);
100 1440000 : auto * node_out_j = _subchannel_mesh.getChannelNode(jj_ch, iz);
101 1440000 : auto rho_i = (*_rho_soln)(node_in_i);
102 1440000 : auto rho_j = (*_rho_soln)(node_in_j);
103 1440000 : auto Si = (*_S_flow_soln)(node_in_i);
104 1440000 : auto Sj = (*_S_flow_soln)(node_in_j);
105 : auto u_star = 0.0;
106 : // figure out donor axial velocity
107 1440000 : if (_Wij(i_gap, iz) > 0.0)
108 383928 : u_star = (*_mdot_soln)(node_out_i) / Si / rho_i;
109 : else
110 1056072 : u_star = (*_mdot_soln)(node_out_j) / Sj / rho_j;
111 :
112 1440000 : crossflow_term +=
113 1440000 : _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz) * u_star;
114 :
115 1440000 : turbulent_term += _WijPrime(i_gap, iz) * (2 * (*_mdot_soln)(node_out) / rho_in / S -
116 1440000 : (*_mdot_soln)(node_out_j) / Sj / rho_j -
117 1440000 : (*_mdot_soln)(node_out_i) / Si / rho_i);
118 1440000 : counter++;
119 : }
120 504000 : turbulent_term *= _CT;
121 :
122 504000 : auto Re = (((*_mdot_soln)(node_in) / S) * Dh_i / mu_in);
123 504000 : auto fi = computeFrictionFactor(Re);
124 504000 : auto ki = k_grid[i_ch][iz - 1];
125 504000 : auto friction_term = (fi * dz / Dh_i + ki) * 0.5 * (std::pow((*_mdot_soln)(node_out), 2.0)) /
126 504000 : (S * (*_rho_soln)(node_out));
127 504000 : auto gravity_term = _g_grav * (*_rho_soln)(node_out)*dz * S;
128 504000 : auto DP = std::pow(S, -1.0) * (time_term + mass_term1 + mass_term2 + crossflow_term +
129 504000 : turbulent_term + friction_term + gravity_term); // Pa
130 504000 : _DP_soln->set(node_out, DP);
131 : }
132 12000 : }
133 1494 : }
134 :
135 : void
136 0 : TriInterWrapper1PhaseProblem::computeh(int iblock)
137 : {
138 0 : unsigned int last_node = (iblock + 1) * _block_size;
139 0 : unsigned int first_node = iblock * _block_size + 1;
140 :
141 0 : if (iblock == 0)
142 : {
143 0 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
144 : {
145 0 : auto * node = _subchannel_mesh.getChannelNode(i_ch, 0);
146 0 : auto h_out = _fp->h_from_p_T((*_P_soln)(node) + _P_out, (*_T_soln)(node));
147 0 : if (h_out < 0)
148 : {
149 0 : mooseError(
150 0 : name(), " : Calculation of negative Enthalpy h_out = : ", h_out, " Axial Level= : ", 0);
151 : }
152 0 : _h_soln->set(node, h_out);
153 : }
154 : }
155 :
156 0 : for (unsigned int iz = first_node; iz < last_node + 1; iz++)
157 : {
158 0 : auto z_grid = _subchannel_mesh.getZGrid();
159 0 : auto dz = z_grid[iz] - z_grid[iz - 1];
160 0 : auto heated_length = _subchannel_mesh.getHeatedLength();
161 0 : auto unheated_length_entry = _subchannel_mesh.getHeatedLengthEntry();
162 :
163 : Real gedge_ave = 0.0;
164 : Real mdot_sum = 0.0;
165 : Real si_sum = 0.0;
166 0 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
167 : {
168 0 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
169 0 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
170 : {
171 0 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
172 0 : auto Si = (*_S_flow_soln)(node_in);
173 0 : auto mdot_in = (*_mdot_soln)(node_in);
174 0 : mdot_sum = mdot_sum + mdot_in;
175 0 : si_sum = si_sum + Si;
176 : }
177 : }
178 0 : gedge_ave = mdot_sum / si_sum;
179 :
180 0 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
181 : {
182 0 : const Real & pitch = _subchannel_mesh.getPitch();
183 0 : const Real & assembly_diameter = _subchannel_mesh.getSideX();
184 0 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
185 0 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, iz);
186 0 : auto mdot_in = (*_mdot_soln)(node_in);
187 0 : auto h_in = (*_h_soln)(node_in); // J/kg
188 0 : auto volume = dz * (*_S_flow_soln)(node_in);
189 0 : auto mdot_out = (*_mdot_soln)(node_out);
190 0 : auto h_out = 0.0;
191 : Real sumWijh = 0.0;
192 : Real sumWijPrimeDhij = 0.0;
193 : Real e_cond = 0.0;
194 :
195 : Real added_enthalpy;
196 0 : if (z_grid[iz] > unheated_length_entry && z_grid[iz] <= unheated_length_entry + heated_length)
197 0 : added_enthalpy = ((*_q_prime_soln)(node_out) + (*_q_prime_soln)(node_in)) * dz / 2.0;
198 : else
199 : added_enthalpy = 0.0;
200 :
201 : // compute the sweep flow enthalpy change
202 0 : auto subch_type = _subchannel_mesh.getSubchannelType(i_ch);
203 : Real sweep_enthalpy = 0.0;
204 :
205 0 : if (subch_type == EChannelType::EDGE || subch_type == EChannelType::CORNER)
206 : {
207 0 : const Real & pitch = _subchannel_mesh.getPitch();
208 0 : const Real & assembly_diameter = _subchannel_mesh.getSideX();
209 : const Real & wire_lead_length = 0.0;
210 : const Real & wire_diameter = 0.0;
211 0 : auto gap = _tri_sch_mesh.getDuctToPinGap();
212 0 : auto w = assembly_diameter + gap;
213 : auto theta =
214 0 : std::acos(wire_lead_length /
215 0 : std::sqrt(std::pow(wire_lead_length, 2) +
216 0 : std::pow(libMesh::pi * (assembly_diameter + wire_diameter), 2)));
217 : // in/out channels for i_ch
218 0 : auto sweep_in = _tri_sch_mesh.getSweepFlowChans(i_ch).first;
219 0 : auto * node_sin = _subchannel_mesh.getChannelNode(sweep_in, iz - 1);
220 0 : auto cs_t = 0.75 * std::pow(wire_lead_length / assembly_diameter, 0.3);
221 0 : auto ar2 = libMesh::pi * (assembly_diameter + wire_diameter) * wire_diameter / 4.0;
222 0 : auto a2p = pitch * (w - assembly_diameter / 2.0) -
223 0 : libMesh::pi * std::pow(assembly_diameter, 2) / 8.0;
224 0 : auto Sij_in = dz * gap;
225 : auto Sij_out = dz * gap;
226 0 : auto wsweep_in = gedge_ave * cs_t * std::pow((ar2 / a2p), 0.5) * std::tan(theta) * Sij_in;
227 0 : auto wsweep_out = gedge_ave * cs_t * std::pow((ar2 / a2p), 0.5) * std::tan(theta) * Sij_out;
228 0 : auto sweep_hin = (*_h_soln)(node_sin);
229 0 : auto sweep_hout = (*_h_soln)(node_in);
230 0 : sweep_enthalpy = (wsweep_in * sweep_hin - wsweep_out * sweep_hout);
231 : }
232 :
233 : // Calculate sum of crossflow into channel i from channels j around i
234 : unsigned int counter = 0;
235 0 : for (auto i_gap : _subchannel_mesh.getChannelGaps(i_ch))
236 : {
237 0 : auto chans = _subchannel_mesh.getGapChannels(i_gap);
238 : unsigned int ii_ch = chans.first;
239 : // i is always the smallest and first index in the mapping
240 : unsigned int jj_ch = chans.second;
241 0 : auto * node_in_i = _subchannel_mesh.getChannelNode(ii_ch, iz - 1);
242 0 : auto * node_in_j = _subchannel_mesh.getChannelNode(jj_ch, iz - 1);
243 : // Define donor enthalpy
244 : auto h_star = 0.0;
245 0 : if (_Wij(i_gap, iz) > 0.0)
246 0 : h_star = (*_h_soln)(node_in_i);
247 0 : else if (_Wij(i_gap, iz) < 0.0)
248 0 : h_star = (*_h_soln)(node_in_j);
249 : // take care of the sign by applying the map, use donor cell
250 0 : sumWijh += _subchannel_mesh.getCrossflowSign(i_ch, counter) * _Wij(i_gap, iz) * h_star;
251 0 : sumWijPrimeDhij += _WijPrime(i_gap, iz) * (2 * (*_h_soln)(node_in) - (*_h_soln)(node_in_j) -
252 0 : (*_h_soln)(node_in_i));
253 0 : counter++;
254 :
255 : // compute the radial heat conduction through gaps
256 0 : auto subch_type_i = _subchannel_mesh.getSubchannelType(ii_ch);
257 0 : auto subch_type_j = _subchannel_mesh.getSubchannelType(jj_ch);
258 : Real dist_ij = pitch;
259 :
260 0 : if (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::EDGE)
261 : {
262 0 : dist_ij = pitch;
263 : }
264 0 : else if ((subch_type_i == EChannelType::CORNER && subch_type_j == EChannelType::EDGE) ||
265 0 : (subch_type_i == EChannelType::EDGE && subch_type_j == EChannelType::CORNER))
266 : {
267 0 : dist_ij = pitch;
268 : }
269 : else
270 : {
271 0 : dist_ij = pitch / std::sqrt(3);
272 : }
273 :
274 0 : auto Sij = dz * _subchannel_mesh.getGapWidth(i_gap);
275 0 : auto thcon_i = _fp->k_from_p_T((*_P_soln)(node_in_i), (*_T_soln)(node_in_i));
276 0 : auto thcon_j = _fp->k_from_p_T((*_P_soln)(node_in_j), (*_T_soln)(node_in_j));
277 : auto shape_factor =
278 0 : 0.66 * (pitch / assembly_diameter) *
279 0 : std::pow((_subchannel_mesh.getGapWidth(i_gap) / assembly_diameter), -0.3);
280 0 : if (ii_ch == i_ch)
281 : {
282 0 : e_cond += 0.5 * (thcon_i + thcon_j) * Sij * shape_factor *
283 0 : ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
284 : }
285 : else
286 : {
287 0 : e_cond += -0.5 * (thcon_i + thcon_j) * Sij * shape_factor *
288 0 : ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) / dist_ij;
289 : }
290 : }
291 :
292 : // compute the axial heat conduction between current and lower axial node
293 0 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz);
294 0 : auto * node_in_j = _subchannel_mesh.getChannelNode(i_ch, iz - 1);
295 0 : auto thcon_i = _fp->k_from_p_T((*_P_soln)(node_in_i), (*_T_soln)(node_in_i));
296 0 : auto thcon_j = _fp->k_from_p_T((*_P_soln)(node_in_j), (*_T_soln)(node_in_j));
297 0 : auto Si = (*_S_flow_soln)(node_in_i);
298 0 : auto dist_ij = z_grid[iz] - z_grid[iz - 1];
299 :
300 0 : e_cond += 0.5 * (thcon_i + thcon_j) * Si * ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) /
301 : dist_ij;
302 :
303 0 : unsigned int nz = _subchannel_mesh.getNumOfAxialCells();
304 : // compute the axial heat conduction between current and upper axial node
305 0 : if (iz < nz)
306 : {
307 :
308 0 : auto * node_in_i = _subchannel_mesh.getChannelNode(i_ch, iz);
309 0 : auto * node_in_j = _subchannel_mesh.getChannelNode(i_ch, iz + 1);
310 0 : auto thcon_i = _fp->k_from_p_T((*_P_soln)(node_in_i), (*_T_soln)(node_in_i));
311 0 : auto thcon_j = _fp->k_from_p_T((*_P_soln)(node_in_j), (*_T_soln)(node_in_j));
312 0 : auto Si = (*_S_flow_soln)(node_in_i);
313 0 : auto dist_ij = z_grid[iz + 1] - z_grid[iz];
314 0 : e_cond += 0.5 * (thcon_i + thcon_j) * Si * ((*_T_soln)(node_in_j) - (*_T_soln)(node_in_i)) /
315 : dist_ij;
316 : }
317 :
318 : // end of radial heat conduction calc.
319 :
320 0 : h_out =
321 0 : (mdot_in * h_in - sumWijh - sumWijPrimeDhij + added_enthalpy + e_cond + sweep_enthalpy +
322 0 : _TR * _rho_soln->old(node_out) * _h_soln->old(node_out) * volume / _dt) /
323 0 : (mdot_out + _TR * (*_rho_soln)(node_out)*volume / _dt);
324 :
325 0 : if (h_out < 0)
326 : {
327 0 : mooseWarning(name(),
328 : " : Calculation of negative Enthalpy h_out = : ",
329 : h_out,
330 : " Axial Level= : ",
331 : iz);
332 : }
333 0 : _h_soln->set(node_out, h_out); // J/kg
334 : }
335 0 : }
336 0 : }
337 :
338 : void
339 30 : TriInterWrapper1PhaseProblem::externalSolve()
340 : {
341 30 : initializeSolution();
342 30 : _console << "Executing subchannel solver\n";
343 30 : auto P_error = 1.0;
344 30 : unsigned int P_it = 0;
345 30 : unsigned int P_it_max = 2 * _n_blocks;
346 30 : if (_n_blocks == 1)
347 : P_it_max = 1;
348 72 : while (P_error > _P_tol && P_it < P_it_max)
349 : {
350 42 : P_it += 1;
351 42 : if (P_it == P_it_max and _n_blocks != 1)
352 : {
353 0 : _console << "Reached maximum number of axial pressure iterations" << std::endl;
354 0 : _converged = false;
355 : }
356 42 : _console << "Solving Outer Iteration : " << P_it << std::endl;
357 42 : auto P_L2norm_old_axial = _P_soln->L2norm();
358 102 : for (unsigned int iblock = 0; iblock < _n_blocks; iblock++)
359 : {
360 60 : int last_level = (iblock + 1) * _block_size;
361 60 : int first_level = iblock * _block_size + 1;
362 60 : auto T_block_error = 1.0;
363 : auto T_it = 0;
364 60 : _console << "Solving Block: " << iblock << " From first level: " << first_level
365 60 : << " to last level: " << last_level << std::endl;
366 :
367 120 : while (T_block_error > _T_tol && T_it < _T_maxit)
368 : {
369 60 : T_it += 1;
370 60 : if (T_it == _T_maxit)
371 : {
372 0 : _console << "Reached maximum number of temperature iterations for block: " << iblock
373 0 : << std::endl;
374 0 : _converged = false;
375 : }
376 60 : auto T_L2norm_old_block = _T_soln->L2norm();
377 :
378 : // computeWij(iblock);
379 60 : computeWijFromSolve(iblock);
380 :
381 60 : if (_compute_power)
382 : {
383 0 : computeh(iblock);
384 :
385 0 : computeT(iblock);
386 : }
387 :
388 60 : if (_compute_density)
389 48 : computeRho(iblock);
390 :
391 60 : if (_compute_viscosity)
392 48 : computeMu(iblock);
393 :
394 : // We must do a global assembly to make sure data is parallel consistent before we do things
395 : // like compute L2 norms
396 60 : _aux->solution().close();
397 :
398 60 : auto T_L2norm_new = _T_soln->L2norm();
399 60 : T_block_error =
400 60 : std::abs((T_L2norm_new - T_L2norm_old_block) / (T_L2norm_old_block + 1E-14));
401 60 : _console << "T_block_error: " << T_block_error << std::endl;
402 : }
403 : }
404 42 : auto P_L2norm_new_axial = _P_soln->L2norm();
405 42 : P_error =
406 42 : std::abs((P_L2norm_new_axial - P_L2norm_old_axial) / (P_L2norm_old_axial + _P_out + 1E-14));
407 42 : _console << "P_error :" << P_error << std::endl;
408 : }
409 : // update old crossflow matrix
410 30 : _Wij_old = _Wij;
411 30 : _console << "Finished executing subchannel solver\n";
412 30 : _aux->solution().close();
413 :
414 : auto power_in = 0.0;
415 : auto power_out = 0.0;
416 1290 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
417 : {
418 1260 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, 0);
419 1260 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, _n_cells);
420 1260 : power_in += (*_mdot_soln)(node_in) * (*_h_soln)(node_in);
421 1260 : power_out += (*_mdot_soln)(node_out) * (*_h_soln)(node_out);
422 : }
423 :
424 : /// TODO: add a verbose print flag
425 30 : auto Total_surface_area = 0.0;
426 30 : auto mass_flow_in = 0.0;
427 30 : auto mass_flow_out = 0.0;
428 1290 : for (unsigned int i_ch = 0; i_ch < _n_channels; i_ch++)
429 : {
430 1260 : auto * node_in = _subchannel_mesh.getChannelNode(i_ch, 0);
431 1260 : auto * node_out = _subchannel_mesh.getChannelNode(i_ch, _n_cells);
432 1260 : Total_surface_area += (*_S_flow_soln)(node_in);
433 1260 : mass_flow_in += (*_mdot_soln)(node_in);
434 1260 : mass_flow_out += (*_mdot_soln)(node_out);
435 : }
436 30 : auto h_bulk_out = power_out / mass_flow_out;
437 30 : auto T_bulk_out = _fp->T_from_p_h(_P_out, h_bulk_out);
438 :
439 30 : _console << " ======================================= " << std::endl;
440 30 : _console << " ======== Subchannel Print Outs ======== " << std::endl;
441 30 : _console << " ======================================= " << std::endl;
442 30 : _console << "Total Surface Area :" << Total_surface_area << " m^2" << std::endl;
443 30 : _console << "Bulk coolant temperature at outlet :" << T_bulk_out << " K" << std::endl;
444 30 : _console << "Power added to coolant is: " << power_out - power_in << " Watt" << std::endl;
445 30 : _console << "Mass in: " << mass_flow_in << " kg/sec" << std::endl;
446 30 : _console << "Mass out: " << mass_flow_out << " kg/sec" << std::endl;
447 30 : _console << " ======================================= " << std::endl;
448 30 : }
|
