Open moosebuild opened 2 years ago
moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @joshuahansel at Aug 12, 2021 09:14AM MST
@snschune Do you happen to still have that writeup of the derivation of pressure drop? If not, that's ok.
moosebuild
commented
2 years ago moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @joshuahansel at Aug 20, 2021 07:18AM MST
@snschune
I did some runs with this test problem:
My input file (Same as yours with minor edits and new PPs):
nelem = 100
friction_factor = 1e6
area = ${fparse 0.5 * 0.353504}
diameter = ${fparse sqrt(4 * area / pi)}
mfr_nominal = 1
p_out_init = 7e6
p_out_nominal = 7e6
T_in = 300
[GlobalParams]
gravity_vector = '0 0 0'
initial_T = ${T_in}
initial_p = ${p_out_init}
initial_vel = 0
initial_vel_x = 0
initial_vel_y = 0
initial_vel_z = 0
closures = none
rdg_slope_reconstruction = minmod
[]
[FluidProperties]
gamma = 1.3066
molar_mass = 2.016e-3
Ru = 8.3145
[h2]
type = IdealGasFluidProperties
gamma = ${gamma}
molar_mass = ${molar_mass}
k = 0.437
mu = 3e-5
[]
[]
[Materials]
[f_wall_mat]
type = ConstantMaterial
property_name = 'f_D'
derivative_vars = 'rhoA rhouA rhoEA'
value = 0.0
block = 'ch_1 ch_3'
[]
[f_wall_mat_2]
type = ConstantMaterial
property_name = 'f_D'
derivative_vars = 'rhoA rhouA rhoEA'
value = ${friction_factor}
block = 'ch_2'
[]
[]
[Components]
[bc_inlet]
type = InletMassFlowRateTemperature1Phase
input = 'ch_1:in'
m_dot = ${mfr_nominal}
T = ${T_in}
[]
[ch_1]
type = FlowChannel1Phase
position = '0 0 0'
orientation = '1 0 0'
length = 1
n_elems = ${nelem}
A = ${area}
D_h = ${diameter}
fp = h2
[]
[junction_1]
type = VolumeJunction1Phase
volume = 1e-5
connections = 'ch_1:out ch_2:in'
position = '1 0 0'
[]
[ch_2]
type = FlowChannel1Phase
position = '1 0 0'
# position = '0 0 0'
orientation = '1 0 0'
# orientation = '0 1 0'
length = 1
n_elems = ${nelem}
A = ${area}
D_h = ${diameter}
fp = h2
[]
[junction_2]
type = VolumeJunction1Phase
volume = 1e-5
connections = 'ch_2:out ch_3:in'
position = '2 0 0'
[]
[ch_3]
type = FlowChannel1Phase
position = '2 0 0'
# position = '1 0.5 0'
orientation = '1 0 0'
# orientation = '-1 0 0'
length = 1
n_elems = ${nelem}
A = ${area}
D_h = ${diameter}
fp = h2
initial_p = ${p_out_init}
# initial_p = ${p_out_nominal}
[]
[bc_outlet]
type = Outlet1Phase
input = 'ch_3:out'
p = ${p_out_nominal}
# p = ${p_out_init}
[]
[]
[Functions]
[mfr_fn]
type = PiecewiseLinear
x = '0 5'
y = '1e-6 ${mfr_nominal}'
[]
[p_out_fn]
type = PiecewiseLinear
x = '0 10'
y = '${p_out_init} ${p_out_nominal}'
[]
[]
[ControlLogic]
[mfr_cntrl]
type = TimeFunctionComponentControl
component = bc_inlet
parameter = m_dot
function = mfr_fn
[]
# [outlet_pressure_cntrl]
# type = TimeFunctionComponentControl
# component = bc_outlet
# parameter = p
# function = p_out_fn
# []
[]
[Preconditioning]
[smp]
type = SMP
full = true
[]
[]
[Postprocessors]
[friction_factor]
type = FunctionValuePostprocessor
function = ${friction_factor}
[]
[p_in_junction]
type = ScalarVariable
variable = junction_1:p
[]
[pressure_drop_junction_minus_p_out]
type = ParsedPostprocessor
pp_names = 'p_in_junction'
function = 'p_in_junction - ${p_out_nominal}'
[]
[p_in_solution]
type = SideAverageValue
variable = p
boundary = 'ch_2:in'
[]
[p_out_solution]
type = SideAverageValue
variable = p
boundary = 'ch_2:out'
[]
[pressure_drop_solution]
type = ParsedPostprocessor
pp_names = 'p_in_solution p_out_solution'
function = 'p_in_solution - p_out_solution'
[]
[momflux_in]
type = FlowJunctionFlux1Phase
junction = junction_1
boundary = 'ch_2:in'
connection_index = 1
equation = momentum
[]
[momflux_out]
type = FlowJunctionFlux1Phase
junction = junction_2
boundary = 'ch_2:out'
connection_index = 0
equation = momentum
[]
[momflux_change]
type = ParsedPostprocessor
pp_names = 'momflux_in momflux_out'
function = '(momflux_in - momflux_out) / ${area}'
[]
[T_in]
type = SideAverageValue
variable = T
boundary = 'ch_2:in'
[]
[T_out]
type = SideAverageValue
variable = T
boundary = 'ch_2:out'
[]
[T_change]
type = ParsedPostprocessor
pp_names = 'T_in T_out'
function = 'T_out - T_in'
[]
[]
[Executioner]
type = Transient
start_time = 0
end_time = 60
# dt = 0.01
[TimeStepper]
type = IterationAdaptiveDT
dt = 1e-5
[]
steady_state_detection = true
steady_state_start_time = 10
n_max_nonlinear_pingpong = 10
timestep_tolerance = 1e-8
petsc_options_iname = '-pc_type'
petsc_options_value = ' lu '
#line_search = l2
nl_rel_tol = 1e-6
nl_abs_tol = 5e-5
nl_max_its = 15
l_tol = 1e-4
l_max_its = 15
automatic_scaling = true
# dtmin = 1e-5
dtmax = 1
error_on_dtmin = false
[]
[Outputs]
[csv]
type = CSV
file_base = 'pressure_drop'
show = 'friction_factor pressure_drop_junction_minus_p_out pressure_drop_solution momflux_change'
execute_on = 'FINAL'
[]
[console]
type = Console
max_rows = 1
[]
[]My plot script:
from math import pi, sqrt
import matplotlib.pyplot as plt
import numpy as np
from thm_utilities import readCSVFile
mfr = 1.0
A = 0.5 * 0.353504
L = 1.0
TL = 300.0
pR = 7e6
Dh = sqrt(4 * A / pi)
Ru = 8.3144598 # universal gas constant
gamma = 1.3066
M = 2.016e-3 # molar mass
R = Ru / M # specific gas constant
cp = gamma * R / (gamma - 1.0)
cv = cp / gamma
def computePressureDropIsothermalNoDynamic(f):
pL = sqrt(pR**2 + f * mfr**2 * R * TL * L / (Dh * A**2))
return pL - pR
log10f_values = [2, 3, 4, 5, 6]
f_values = [10**log10f for log10f in log10f_values]
dp_analytic = [computePressureDropIsothermalNoDynamic(f) for f in f_values]
data_10elem = readCSVFile("pressure_drop_10elem.csv")
data_100elem = readCSVFile("pressure_drop_100elem.csv")
data_10elem_minmod = readCSVFile("pressure_drop_10elem_minmod.csv")
data_100elem_minmod = readCSVFile("pressure_drop_100elem_minmod.csv")
data_10elem_full = readCSVFile("pressure_drop_10elem_full.csv")
data_100elem_full = readCSVFile("pressure_drop_100elem_full.csv")
plt.figure(figsize=(8, 6))
plt.rc('text', usetex=True)
plt.rc('font', family='sans-serif')
ax = plt.subplot(1, 1, 1)
ax.get_yaxis().get_major_formatter().set_useOffset(False)
plt.xlabel("Friction Factor")
plt.ylabel("Pressure Drop [Pa]")
plt.loglog(f_values, dp_analytic, '-', color='black', marker='', label='Analytic')
plt.loglog(data_10elem['friction_factor'], data_10elem['pressure_drop_junction_minus_p_out'], '-', color='indianred', marker='', label='$p_J - p_{out}$, 10 elements')
plt.loglog(data_10elem['friction_factor'], data_10elem['pressure_drop_solution'], '--', color='indianred', marker='', label='$p_1 - p_N$, 10 elements')
plt.loglog(data_10elem['friction_factor'], data_10elem['momflux_change'] / A, ':', color='indianred', marker='', label='$(\\rho u^2 + p)_L - (\\rho u^2 + p)_R$, 10 elements')
plt.loglog(data_100elem['friction_factor'], data_100elem['pressure_drop_junction_minus_p_out'], '-', color='orange', marker='', label='$p_J - p_{out}$, 100 elements')
plt.loglog(data_100elem['friction_factor'], data_100elem['pressure_drop_solution'], '--', color='orange', marker='', label='$p_1 - p_N$, 100 elements')
plt.loglog(data_100elem['friction_factor'], data_100elem['momflux_change'], ':', color='orange', marker='', label='$(\\rho u^2 + p)_L - (\\rho u^2 + p)_R$, 100 elements')
plt.loglog(data_10elem_minmod['friction_factor'], data_10elem_minmod['pressure_drop_junction_minus_p_out'], '-', color='lightgreen', marker='', label='$p_J - p_{out}$, 10 elements with minmod slope reconstruction')
plt.loglog(data_100elem_minmod['friction_factor'], data_100elem_minmod['pressure_drop_junction_minus_p_out'], '-', color='turquoise', marker='', label='$p_J - p_{out}$, 100 elements with minmod slope reconstruction')
plt.loglog(data_10elem_full['friction_factor'], data_10elem_full['pressure_drop_junction_minus_p_out'], '-', color='cornflowerblue', marker='', label='$p_J - p_{out}$, 10 elements with full slope reconstruction')
plt.loglog(data_100elem_full['friction_factor'], data_100elem_full['pressure_drop_junction_minus_p_out'], '-', color='mediumorchid', marker='', label='$p_J - p_{out}$, 100 elements with full slope reconstruction')
ax.legend(frameon=False, prop={'size':12}, loc='upper left')
plt.tight_layout()
plt.savefig('dp_vs_f.png', dpi=300)Here's a summary of the cases I ran:
p_J - p_out: This is what you did (subtracted provided outlet pressure from the upstream junction pressure)p_1 - p_N: Subtracted last element pressure from first element pressure(rhou^2 + p)_L - (rhou^2 + p)_R: Subtracted downstream momentum flux from upstream momentum flux.I also checked the isothermal assumption that you used in the derivation and found it to be good; the temperature change was O(-3) K. The (rhou^2 + p)_L - (rhou^2 + p)_R sets were used to check the implicit assumption d/dx(rhou^2) = 0, which again, is good. So the reference solution should be essentially correct.
Here's what I get from the results:
Some notes: I think the default configuration that we should suggest to a user is slope reconstruction with limitation (in contrast to our current default, Godunov), like minmod. This problem happened to be well behaved enough that full slope reconstruction could be used without risk. Maybe if we find that a large enough class of problems performs similarly well, full reconstruction should be the default. We currently have Godunov as the default because we didn't have the slope Jacobians yet and at least with slope limiters like minmod, the nonlinear convergence can be significantly deteriorated, at least for problems where it's activated. It's been reported that this can be the case, even when the Jacobians are correctly computed, just due to the high degree of nonlinearity.
Anyway, I still plan to look at all-Mach modifications to HLLC, as I think it has the potential to severely reduce the artificial dissipation seen by the Godunov scheme in the low Mach regime, and probably make noticeable improvements even when using slope reconstruction.
@lindsayad @andrsd @charlise
moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @lindsayad at Aug 20, 2021 08:27AM MST
Does "full" reconstruction correspond to linear interpolation from cell centroid values to the face? If not, what is it?
moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @snschune at Aug 20, 2021 08:50AM MST
Thanks @joshuahansel that is great progress. The other issue we always had was that convergence deteriorated when increasing the friction factor. Do you observe that as well? I'll use your recommendation in my NTP model later today or next week to see what's going on.
moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @joshuahansel at Aug 20, 2021 09:09AM MST
Does "full" reconstruction correspond to linear interpolation from cell centroid values to the face?
Yes, we construct slopes for {p,u,T} using cell-centroid values for those quantities and then interpolate to faces.
The other issue we always had was that convergence deteriorated when increasing the friction factor. Do you observe that as well?
Yes, but I haven't investigated it yet. I'll try it again after we get THM converted to AD, which I expect to be within a week, just to make sure Jacobians aren't to blame. If poor convergence holds (likely), I'm not sure what can be done - do we have reason to believe that convergence shouldn't be poor for such high friction factors? In other words, are we sure that the physics aren't just inherently very nonlinear in this case? When would these friction factors be used? Is this just to model a turbine? If so, maybe we wouldn't need good performance here if we had a functional 0-D turbine model that could give the desired pressure drop (Lise is working on this).
moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @lindsayad at Aug 20, 2021 09:14AM MST
Yes
I'm plainly shocked that this is stable then. Is this only because implicit Euler is being used as the time integrator? What happens if we change the time integrator or decrease the time step? I would think oscillations should appear.
moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @joshuahansel at Aug 20, 2021 09:39AM MST
decrease the time step
Are you saying it should be less stable for smaller time steps? The runs use an adaptive time step size, and it basically ends up being the largest time step size such that it doesn't fail the time step.
change the time integrator
We usually recommend BDF2, unless we're doing some shock problem, in which case we'd be fine with explicit methods. Do you just want me to compare BE and BDF2? I can look at the pressure transients to see if oscillations appear in either.
moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @lindsayad at Aug 20, 2021 09:57AM MST
Yea I'm just used to seeing oscillations appear any time you have a central-differencing-esque discretization for the advection term in a transient simulation. For instance in scalar advection with linear interpolation to faces you would not have any on-diagonal from the advection term and it's notoriously unstable
moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @joshuahansel at Aug 20, 2021 10:04AM MST
Just to be clear, HLLC is still used to compute the fluxes in this case. It's just that the input states are computed from the linear interpolation.
moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @lindsayad at Aug 20, 2021 10:45AM MST
If you are not limiting and you are linearly interpolating, then the input states, e.g. the states on the + and - side of the face, are the exact same right?
But as you say... you are computing the fluxes with HLLC ... So perhaps that helps the stability? E.g. perhaps with the HLLC computation of fluxes the matrix corresponding to advection discretization does still have on-diagonals even with linear interpolation?
moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @joshuahansel at Aug 20, 2021 12:05PM MST
If you are not limiting and you are linearly interpolating, then the input states, e.g. the states on the + and - side of the face, are the exact same right?
Nope, I didn't explain it well. A single slope value is computed for each cell; for cell i, slope_i depends on {u_{i-1}, u_i, u_{i+1}}. On the left side of i+1/2, you use u_i + slope_i * (x_{i+1/2} - x_i), and on the right side you use u_{i+1} + slope_{i+1} * (x_{i+1/2} - x_{i+1}).
moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @lindsayad at Aug 20, 2021 01:03PM MST
This is very interesting. The equations in Moukalled and the limiter implementation in OpenFOAM essentially limit an interpolation between the current cell C value and the "downwind" cell D value (when doing directional interpolations like we do for HLLC we go both ways, so "downwind" is not always downwind). However, what you are doing is limiting the cell center gradient and then using that limited cell center gradient to perform extrapolation. As a result of this when doing full (unlimited) slope reconstruction, your + side face value will indeed be different from your - side face value and both those values will have a larger stencil (e.g. your + side value (I'll associate the + side with cell C) will depend on U, C, D values and your - side values will depend on C, D, and DD values) whereas our central difference interpolation will result in the same + and - values and both stencils will be small, both depending only on C and D values.
I'm curious what reference you guys used for implementing your limiters?
moosebuild
commented
2 years ago Imported from idaholab/thm#214 by @joshuahansel at Aug 20, 2021 01:21PM MST
I'm curious what reference you guys used for implementing your limiters?
I don't have my textbooks unpacked yet, so I can't verify, but my first sources would have been Toro's book that you're familiar with (the MUSCL-Hancock section maybe?) and then LeVeque's book "Finite Volume Methods for Hyperbolic Problem". The latter is the one we reference in our theory manual, but I'm not sure if it does the same thing we do or if it was a generic reference on the subject. Yidong was the originator of this slope reconstruction implementation.
lindsayad
commented
2 years ago @snschune @joshuahansel should we regard this issue as closed? @joshuahansel have you changed the default config to be slope reconstruction with limitation or is the default still Godunov (piecewise constant)?
joshuahansel
commented
2 years ago I have not changed the default in THM to use slope reconstruction. I still need to do this.
lindsayad
commented
2 years ago @grmnptr, see https://github.com/idaholab/moose/issues/19894#issuecomment-999937664. I believe that we have had some discussion about this before, e.g. the difference in limiting surface interpolations vs. limiting gradients. Don't they have distinct schemes for these two concepts in OpenFOAM?
lindsayad
commented
2 years ago @grmnptr what do you think about these two different options for computing face values? E.g. limited surface interpolation vs. computation of "sided" face values using a cell center value and limited cell center gradient? The latter is perhaps only appropriate with an advection scheme that supports having a discontinuous jump in face value? (E.g. KT or HLLC)
grmnptr
commented
2 years ago I think we can explore the second one as well at some point. In some cases we need gradient limiting in the first option as well. For example the limiters we have right now use the gradient in the computation of the r parameter that goes into the limiter. Depending on the limiter we might have to limit that gradient. Yes, in OF they have both gradient limiters and surface interpolation limiters. I think adding the second option should be interesting, especially if we are aiming for compressible fluids/regimes. I am also not 100% certain how those schemes translate to problems with unstructured meshes.
Imported from idaholab/thm#214 by @snschune at Aug 12, 2021 08:47AM MST
Bug Description
@joshuahansel
Pressure drop vs. friction factor for a single channel does not match hand calculations (see excel sheet).
relap-7-dp.xlsx
Pressure is taken from the junction but visual inspection of the pressure trace over the three channels in the problem indicates that the problem is also encountered with the end point pressures.
Bug Impact
Potentially incorrect results when distributed pressure drops are used.
Steps to Reproduce Bug