You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
where $\beta$ is the turbulent mixing parameter or thermal transfer coefficient and $\bar{G}$ is the average mass flux of the adjacent subchannels. The $\beta$ term is the tuning parameter for the mixing model. Physically, it is a non-dimensional coefficient that represents the ratio of the lateral mass flux due to mixing to the axial mass flux. It is used to model the effect of the unresolved scales of motion that are produced through the averaging process. In single-phase flow no net mass exchange occurs, both momentum and energy are exchanged between subchannels, and their rates of exchange are characterized in terms of hypothetical turbulent interchange flow rates ($w_{ij}^{'H},w_{ij}^{'M}$) [!cite](TODREAS), for enthalpy and momentum respectively. The approximation that the rate of turbulent exchange for energy and momentum are equal is adopted: $w'_{ij} = w_{ij}^{'H} = w_{ij}^{'M}$.
121
+
where $\beta$ is the turbulent mixing parameter or thermal transfer coefficient and $\bar{G}$ is the average mass flux of the adjacent subchannels. The $\beta$ term is the tuning parameter for the mixing model. Physically, it is a non-dimensional coefficient that represents the ratio of the lateral mass flux due to mixing to the axial mass flux. It is used to model the effect of the unresolved scales of motion that are produced through the averaging process. In single-phase flow no net mass exchange occurs, both momentum and energy are exchanged between subchannels, and their rates of exchange are characterized in terms of hypothetical turbulent interchange flow rates ($w_{ij}^{'H},w_{ij}^{'M}$) [!cite](TODREAS), for enthalpy and momentum respectively. The approximation that the rate of turbulent exchange for energy and momentum are related as follows is adopted: $w'_{ij} = w_{ij}^{'H} = w_{ij}^{'M} / C_T$.
122
122
123
123
Additional turbulent mixing parameters are implemented as follows:
124
124
@@ -136,7 +136,6 @@ After calibrating the turbulent diffusion coefficient $\beta$ we turned our atte
136
136
137
137
For quadrilateral assemblies: $C_{T} = 2.6$, $\beta = 0.006$ [!cite](kyriakopoulos2022development).
138
138
139
-
140
139
## Discretization
141
140
142
141
The collocated discretization of the variables is presented in [fig:dis] . $i,j$ are the subchannel indexes. $ij$ is the name of the gap between subchannels $i,j$. $k$ is the index in the axial direction.
@@ -269,11 +268,11 @@ This is the default algorithm, where the unknown flow variables are calculated i
269
268
270
269
#### Implicit segregated
271
270
272
-
In this case, the governing equations are recast in matrix form and the flow variables are calculated by solving the corresponding system. This means that variables are retrieved concurrently for the whole block. Otherwise, the solution algorithm is the same as in the default method.
271
+
In this case, the governing mass, axial momentum and crossflow momentum, equations are recast in matrix form and the flow variables are calculated by solving the corresponding system. This means that variables are retrieved concurrently for the whole block. Otherwise, the solution algorithm is the same as in the default explicit method.
273
272
274
273
#### Implicit
275
274
276
-
In this case, the conservation equations are recast in matrix form and combined into a single system. The user can decide whether or not they will include the enthalpy calculation in the matrix formulation. The flow variables are calculated by solving that big system to retrieve all the unknows at the same time instead of one by one, and on all the nodes of the block: $\vec{\dot{m}}, \vec{P}, \vec{w_{ij}}, \vec{h}$. The solution algorithm is the same as in the default method, but the solver used in this version is a fixed point iteration instead of a Newton method. The system looks like this:
275
+
In this case, the governing mass, axial momentum and crossflow momentum conservation equations are recast in matrix form and combined into a single system. The system of all the subchannel equations looks like this:
277
276
278
277
\begin{equation}
279
278
\begin{bmatrix}
@@ -297,5 +296,93 @@ In this case, the conservation equations are recast in matrix form and combined
297
296
\end{bmatrix}
298
297
\end{equation}
299
298
300
-
As soon as the big matrix is constructed, the solver will calculate cross-flow resistances to achieve realizable solutions. This is done by an initial calculation of the axial mass flows, crossflows, and the sum of crossflows.
301
-
The defined relaxation factor is the average crossflow divided by the maximum and added to 0.5. The added cross-flow resistance is a blending of a current value and the previous one using the relaxation factor calculated above. The current value is the maximum of the sum of cross-flows per subchannel over the minimum axial mass-flow rate. The added cross-flow resistance is added to the diagonal of $M_{ww}$ matrix.
299
+
Since the enthalpy governing equations are uncoupled from the other equations in this otherwise monolithic system (enthalpy is coupled to the flow equations via the fluid properties update), it makes sense to lag the enthalpy solution and solve for it separately. The flow variables are calculated by solving that big system (without the enthalpy) to retrieve all the unknowns at the same time instead of one by one, and on all the nodes of the block: $\vec{\dot{m}}, \vec{P}, \vec{w_{ij}}$. The solution algorithm is the same as in the default method and the solver used is PETSc KSPSolve.
300
+
301
+
As soon as the big matrix is constructed, the solver will calculate cross-flow resistances to maintain realizability. A distinctive feature of this method is the introduction of a *weak relaxation* logic that stabilizes and accelerates convergence of the coupled $mass flow: (\dot{\mathbf{m}})$, $pressure: (\mathbf{P})$, and $crossflow:(\mathbf{w}_{ij})$ fields in a $Q{=}3$ block-nested linear system with matrix blocks $M_{ij}$ and right-hand-side blocks $\mathbf{b}_i$ that represent the individual governing equations. Note that the solution is influenced by the stabilization method and its coefficients.
302
+
303
+
#### 1. Fast scale estimates
304
+
305
+
From the axial- and cross-momentum rows, the code forms quick, diagonally preconditioned estimates:
with small safeguards $\varepsilon_p,\varepsilon_W\sim 10^{-10}$ to avoid division by zero. Using $\hat{\mathbf W}$, the per-channel crossflow sum $\sum_{j} w_{ij}$ is assembled into a vector $\mathrm{sumw_{ij}}_{\mathrm{loc}}$.
r = \frac{r_{\mathrm{base}}}{\max(S_{\max}, \varepsilon)} + 0.5,\qquad \varepsilon\sim10^{-10}.
333
+
\end{equation}
334
+
The +0.5 offset biases toward mild under-relaxation.
335
+
336
+
#### 3. Crossflow resistance inflation
337
+
338
+
A cross-coupling resistance is estimated and smoothed:
339
+
\begin{equation}
340
+
\begin{aligned}
341
+
\tilde K &= \frac{S_{\max}}{m_{\min}}, &
342
+
K^\star &= 0.9\,\tilde K + 0.1\,K_{\text{old}}, &
343
+
K &= r\,K^\star.
344
+
\end{aligned}
345
+
\end{equation}
346
+
After smoothing, the provisional crossflow resistance $K$ is mapped through a piecewise lower-bound function that enforces minimum safe damping levels in specific ranges.
347
+
348
+
\begin{equation}
349
+
K \rightarrow
350
+
\begin{cases}
351
+
K , & K >= 10, \\
352
+
1.0, & 1 \leq K < 10, \\
353
+
0.5, & 0.1 \leq K < 1, \\
354
+
\frac{1}{3}, & 0.01 \leq K < 0.1, \\
355
+
0.1, & 0.001 \leq K < 0.01, \\
356
+
K, & K < 10^{-3}.
357
+
\end{cases}
358
+
\end{equation}
359
+
360
+
This mapping acts as a {snap-up} rule for the crossflow resistance $K$ over the range $[10^{-3}, 10]$:
361
+
it raises $K$ out of weak-damping intervals but leaves very small and very large
362
+
values unchanged. The purpose is to maintain numerical stability and adequate
363
+
diagonal dominance in the cross-momentum equations without introducing full quantization or "bucketing".
364
+
365
+
Finally, $K$ is added to the diagonal of the cross-momentum block,
366
+
\begin{equation}
367
+
M_{ww} \;\leftarrow\; M_{ww} + K\,I,
368
+
\end{equation}
369
+
thereby increasing diagonal dominance and improving conditioning for the crossflow equations. Note that this treatment does influence the cross-flow distribution solution.
370
+
371
+
#### 4. Per-equation under-relaxation
372
+
373
+
Classical linear under-relaxation is applied automatically and separately to each equation $f\in\{\mathbf m,\mathbf p,\mathbf W\}$ using factors
This standard construction ensures that solving the modified linear system yields the *under-relaxed* update for equation $f$. In practice, only $\mathbf W$ is strongly damped, while $\mathbf m$ and $\mathbf p$ can be solved without additional damping. This relaxation happens inside the temperature loop.
385
+
386
+
#### 5. Net effect
387
+
388
+
The combination of (i) safeguarded scale estimation, (ii) adaptive, time smoothed, and piecewise snapped added crossflow resistance, and (iii) selective under-relaxation produces a more diagonally dominant and robust nested solve that tolerates rapid changes in crossflow while preserving good convergence properties for mass flow and pressure.
0 commit comments