Add option to advect enthalpy instead of cp T for non constant cp fluids #29569

GiudGiud · 2024-12-18T18:59:48Z

We already had a correction in the time derivative for non-constant cp. But @freiler spotted that we did not in the advection term.

This should eventually be the default. I ll do that when I transition the nonlinear fluid heat transfer Physics to use this formulation

tanoret · 2024-12-18T19:05:04Z

Please add realistic testing or more extended verification against analytical solutions. Let's try avoiding the same issues that we had in the past. Thank you!

tanoret · 2024-12-18T19:26:11Z

Please add the same tests than for the linear problem:

1D problem with cp=a+bT and verification against analytical solution
2D problem with wall heating and verifying enthalpy balances (h_dot_out = h_dot_in + q''_wall * A_wall)
Equivalent problem in 3D and test against the previous 2D problem (this is important to see that we are not blowing up the size of the AD container by advecting enthalpy)
Body with variable cp and rho cooling at a fixed heat rate and verify against the analytical solution (this is verifying the implementation of the time derivative)

Thank you!

GiudGiud · 2024-12-18T19:37:16Z

we should be good for this one

Equivalent problem in 3D and test against the previous 2D problem (this is important to see that we are not blowing up the size of the AD container by advecting enthalpy)

whether rho cp T or h, you have the same dependencies on pressure and temperature. We are not adding new ones by switching to h

For the time derivative we probably want to do more work than this if we want a conservative scheme.
I ll see what I can do for 1D / 2D

freiler · 2024-12-18T20:00:23Z

We can use the 1D and 2D cases I have for enthalpy to eventually check the accuracy of the new kernel. In addition, I will also build a transient case to test the time derivative kernel.

GiudGiud · 2024-12-18T21:38:34Z

@freiler do you want to test this PR to provide 1 and 2? Or should I use your inputs?

freiler · 2024-12-18T21:40:32Z

I can push the input for the 1D and 2D cases to your branch so we at least run the same test cases and compare.

moosebuild · 2024-12-19T00:16:11Z

Job Coverage, step Generate coverage on b22431f wanted to post the following:

Framework coverage

Coverage did not change

Modules coverage

Fluid properties

	e193a3	#29569 b22431
	Total	Total	+/-	New
Rate	85.11%	85.13%	+0.02%	100.00%
Hits	7850	7863	+13	14
Misses	1373	1373	-	0

Diff coverage report

Full coverage report

Navier stokes

	e193a3	#29569 b22431
	Total	Total	+/-	New
Rate	84.72%	84.74%	+0.02%	95.83%
Hits	17494	17527	+33	46
Misses	3156	3157	+1	2

Diff coverage report

Full coverage report

Full coverage reports

Reports

This comment will be updated on new commits.

moosebuild · 2024-12-19T01:03:27Z

Job Documentation, step Docs: sync website on b22431f wanted to post the following:

View the site here

This comment will be updated on new commits.

refs idaholab#12174

- balance to 1e-8 - error to T_out in analytical solution down to 5e-6

GiudGiud · 2024-12-24T22:25:20Z

I think the transient case should be done in a PR on time integration. As explained in the thread on slack, I dont expect the time derivative discretization rho_new cp_new T_new - rho_new cp_new T_old to be conservative.

I also do not think we would be able to distinguish this time discretization conservation error from the regular order error of a first order time integrator

lindsayad · 2025-01-07T20:09:13Z

@freiler @tanoret can one of you review this please?

grmnptr

Looks good to me! - did the code review. I think the test cases are also convincing.

INSFVEnergyTimeDerivative should be correct if you pass the new functor to it.
WCNSFVEnergyTimeDerivative should also work I think.

Could you please add a test to verify that they work?

grmnptr · 2025-01-11T19:59:37Z

modules/navier_stokes/src/functormaterials/INSFVEnthalpyFunctorMaterial.C

@@ -23,7 +24,7 @@ INSFVEnthalpyFunctorMaterial::validParams()
  params.addClassDescription(
      "This is the material class used to compute enthalpy for the "
      "incompressible/weakly-compressible finite-volume implementation of the Navier-Stokes "
-      "equations. Note that this class assumes that cp is a constant");
+      "equations.");


With these changes this will look similar to @freiler 's functor material. Do you see value in having only one? Or should we keep two? I think if it can be done cleanly, we should have only one.

I havent been keeping track with his PR but I think his functor material might be focusing on computing T from h rather than h from T

If the class assumes that cp is constant, please add it to the name of the functormaterial, e.g., INSFVConstantCpEnthalpyFunctorMaterial or something. Otherwise, this is can lead to errors

the assumptions depend on the parameter

grmnptr · 2025-01-11T20:04:36Z

modules/navier_stokes/src/functormaterials/INSFVEnthalpyFunctorMaterial.C

-  addFunctorProperty<ADReal>(
-      "rho_cp_temp", [&rho_h](const auto & r, const auto & t) -> ADReal { return rho_h(r, t); });
+    addFunctorProperty<ADReal>(
+        "rho_cp_temp", [&rho_h](const auto & r, const auto & t) -> ADReal { return rho_h(r, t); });


If we are making this more general, I'd not call it rho_cp_temp anymore. Most of the code fetches this under the hood so I hope not a lot of input files need changing.

ok let's just not define these functors anymore.

added the transient test too

- stop defining cpT and rho cp T functors - add a transient test + add an exception test

tanoret

Thank you for the work. I am suggesting a couple of changes in the tests. Please let me know if further clarification is needed for anything. Thanks again!

tanoret · 2025-01-12T01:16:40Z

modules/fluid_properties/src/fluidproperties/NaClFluidProperties.C

+    Real pressure, Real temperature, Real & cp, Real & dcp_dp, Real & dcp_dT) const
+{
+  // Correlation needs pressure in bar
+  Real pbar = pressure * 10.0e-5;


Conversion to bar should be 1e-5 (1 bar = 1e5 Pa)

very true. this is like this in the original property :O ( I copy pasted that line)

tanoret · 2025-01-12T01:17:46Z

modules/fluid_properties/src/fluidproperties/NaClFluidProperties.C

+  // Halite isobaric heat capacity
+  cp = 1148.81 + 0.551548 * (Tc - Tt) + 2.64309e-4 * (Tc - Tt) * (Tc - Tt) + r3 * pbar +
+       r4 * pbar * pbar;
+  dcp_dp = r3 * 10.e-5 + 2 * r4 * pbar * 10.e-5;


Conversion to bar may also be wrong here

tanoret · 2025-01-12T01:19:06Z

modules/fluid_properties/src/fluidproperties/NaClFluidProperties.C

  return 1148.81 + 0.551548 * (Tc - Tt) + 2.64309e-4 * (Tc - Tt) * (Tc - Tt) + r3 * pbar +
         r4 * pbar * pbar;
 }

+void
+NaClFluidProperties::cp_from_p_T(
+    Real pressure, Real temperature, Real & cp, Real & dcp_dp, Real & dcp_dT) const


Could we please add unit test to compare the computed cp against tabulated values for a few points in the input sapce?

they already exist, the dependence on pressure is just so small it did not catch it

modules/navier_stokes/include/functormaterials/INSFVEnthalpyFunctorMaterial.h

modules/navier_stokes/src/functormaterials/INSFVEnthalpyFunctorMaterial.C

modules/navier_stokes/test/tests/finite_volume/wcns/materials/1d_test_cpT.i

tanoret · 2025-01-12T01:26:08Z

modules/navier_stokes/test/tests/finite_volume/wcns/materials/1d_test_cpT.i

+  []
+  [T_analytical]
+    type = Receiver
+    default = ${fparse (-A_cp+sqrt(A_cp^2-2*B_cp*(-q_source/rho/bulk_u*L-A_cp*T_in-B_cp/2*T_in*T_in)))/B_cp}


I verified by hand calculation that this analytical solution is correct. However, please add a couple comments of the steps of the derivation for future reference.

Ramiro did the derivation. @freiler do you want to add that as a suggestion?

tanoret · 2025-01-12T01:26:34Z

modules/navier_stokes/test/tests/finite_volume/wcns/materials/1d_test_cpT.i

@@ -0,0 +1,231 @@
+L = 30


This is an excellent test! Exactly the kind of things we need. Thank you!

tanoret · 2025-01-12T01:28:28Z

modules/navier_stokes/test/tests/finite_volume/wcns/materials/enthalpy_computation.i

@@ -0,0 +1,101 @@
+[Mesh]


What is the purpose of this test? Let's compute the enthalpies by hand and compare the code result against the hand calculation. That way we know that things are being done properly

it's the original regression test. I ll add the hand calc

done. Had to change to PbBiFPs because Water97 is quite complex in terms of computing a specific fluid property

tanoret · 2025-01-12T01:30:12Z

modules/navier_stokes/test/tests/finite_volume/wcns/materials/2d-steady-wall-balance.i

+  [balance_in_percent]
+    type = ParsedPostprocessor
+    expression = '(H_out + H_in - Q - Q_wall) / H_in * 100'
+    pp_names = 'H_in H_out Q Q_wall'


This test is currently formulated as a non-regression test. However, this is the only postprocessor that we actually need to check. Let's please compare in the test this pp against 0 with a certain tolerance

comparing stuff to 0 (or to this small value) makes it awkward with the tolerances (relative does not work). I'd almost want to compare the 4 (far from zero) terms without comparing that one? (though it works right now so no need to change it)

modules/navier_stokes/test/tests/finite_volume/wcns/materials/1d_test_cpT.i

modules/navier_stokes/test/tests/finite_volume/wcns/materials/2d-steady-wall-balance.i

modules/navier_stokes/test/tests/finite_volume/wcns/materials/tests

modules/navier_stokes/test/tests/finite_volume/wcns/materials/gold/1d_test_cpT_out.csv

- rename a boolean - rename a postprocessor - fix conversion to bar in NaCl cp computation - simplify enthalpy computation test to add a "hand" calc of the enthalpy Co-authored-by: Mauricio Tano <mauriciotano@gmail.com>

GiudGiud self-assigned this Dec 18, 2024

GiudGiud force-pushed the PR_non_constant_cp branch from bb2c1a6 to ac93fc1 Compare December 18, 2024 19:37

GiudGiud force-pushed the PR_non_constant_cp branch from ac93fc1 to c488692 Compare December 18, 2024 21:37

GiudGiud force-pushed the PR_non_constant_cp branch from 0a0aa7d to f0a01b1 Compare December 19, 2024 21:42

GiudGiud added 2 commits December 19, 2024 22:42

Add an option to use a FP UO to compute the enthalpy

c4dbe49

refs idaholab#12174

Add a test for non-constant cp

02dc1c9

GiudGiud force-pushed the PR_non_constant_cp branch from f0a01b1 to 97db47a Compare December 19, 2024 21:44

idaholab deleted a comment from moosebuild Dec 19, 2024

GiudGiud force-pushed the PR_non_constant_cp branch from 97db47a to d0d9f2f Compare December 19, 2024 22:51

Add an option to input h directly so we can use user-input FPs

4c3f05c

GiudGiud force-pushed the PR_non_constant_cp branch from d0d9f2f to 4d8a100 Compare December 20, 2024 01:40

GiudGiud added 2 commits December 20, 2024 13:44

Add cp derivatives to NaCl, was done for testing purposes

fb7cf27

Add 1D case showing:

5b7a658

- balance to 1e-8 - error to T_out in analytical solution down to 5e-6

GiudGiud force-pushed the PR_non_constant_cp branch 2 times, most recently from 6879736 to 5b7a658 Compare December 20, 2024 15:30

moosebuild added the PR: Failed but allowed label Dec 20, 2024

idaholab deleted a comment from moosebuild Dec 24, 2024

GiudGiud added the PR: Ready for review/merge label Dec 24, 2024

moosebuild removed the PR: Failed but allowed label Dec 24, 2024

GiudGiud marked this pull request as ready for review December 24, 2024 22:25

GiudGiud requested review from lindsayad, grmnptr and joshuahansel as code owners December 24, 2024 22:25

Add a 2D case showing the energy balance in 2D with wall heating

367bbf8

GiudGiud force-pushed the PR_non_constant_cp branch from 76e7282 to 367bbf8 Compare December 26, 2024 20:48

grmnptr reviewed Jan 11, 2025

View reviewed changes

Address Peter's review:

c9faf9e

- stop defining cpT and rho cp T functors - add a transient test + add an exception test

GiudGiud force-pushed the PR_non_constant_cp branch from dfcac83 to c9faf9e Compare January 12, 2025 01:12

tanoret suggested changes Jan 12, 2025

View reviewed changes