Debug: Database connection successful
You are not logged in.
fvModels is a control file that provides "volume" heating of our 240 meter heating system for the Optical Plane vessel.
Here is what it looks like now. We will likely have to make changes to allow it to run with the new mesh...
cat constant/fvModels
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / Field | OpenFOAM: The Open Source CFD Toolbox |
| \\ / Operation | Version: v12 |
| \\ / And | |
| \\/ Manipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object fvModels;
}uniformIntakeHeater
{
type semiImplicitSource;
// type bananaDoesNotExist;
active yes;select cellSet;
cellSet heaterCells;volumeMode absolute; // "absolute" = total quantity over the selected region
sources
{
e
{
explicit
{
type constant;
value 4e7;
}
implicit 0; // 40 MW total over 24,000 cells
}
}
}// Adding new JSON to implement washer to restrict throat
washerRestriction
{
type porosityForce;
active yes;porosityForceCoeffs
{
// IMPORTANT: porosityForce expects selectionMode (cellZone is typical)
selectionMode cellZone;
cellZone washerZone;type DarcyForchheimer;
DarcyForchheimerCoeffs
{
// Start “gentle”, then ramp up by decades until leakage is negligible.
// Strong resistance in x, minimal in y/z: Start 1e8 increment to 1e12
d (1e9 0 0);
f (0 0 0);// Align resistance with x-axis
coordinateSystem
{
type cartesian;
origin (0 0 0);
coordinateRotation
{
type axesRotation;
e1 (1 0 0);
e2 (0 1 0);
}
}
}
}
}
Right off the bat, I can see that we need to remove that washer section.
While we're waiting for ChatGPT5.2 to get back to me, I can toss in an explanation for the banana...
bananaDoesNotExist;
We used this to make sure the fluid solver was picking up fvModels correctly. The confirmation occurred in the form of an error message reporting that bananaDoesNotExist was not found. We saw this error message when we started the case and it hiccuped due to the incorrect name.
Looking ahead, I'm guessing we are going to have to implement cellSet in our new mesh.
That will require adding another control file, and that may be in the system folder which we have not yet addressed, other than for controlDict.
Here is ChatGPT5.2's evaluation: ( I've omitted the confirmation that the washer section needs to go. That section is removed)
One more critical check (because your new mesh is brand new)
Your heater model targets:
select cellSet;
cellSet heaterCells;That means the case will only work if heaterCells exists in the new mesh.
In your new foamWorkroom mesh created by blockMesh, you currently have only:
constant/polyMesh/{boundary,faces,neighbour,owner,points}
No sets/ directory was listed, so you almost certainly do NOT yet have heaterCells. << sets was added
How to confirm
From the case directory:
ls constant/polyMesh/sets
If that directory doesn’t exist, or there’s no heaterCells, then before running you must create it (typically via topoSet using system/topoSetDict, or by switching the model to “all cells” selection).
The two clean ways forward
Option A (recommended, preserves your proven config):
Create heaterCells for the new pipe mesh using topoSet (and any other sets you rely on, like heatedWallFaces, etc.).
This keeps your fvModels exactly compatible with what you already debugged.
The next step is to look at topoSetDict in system on the running case.
(th)
Offline
Like button can go here
ChatGPT5.2 recommended a small adjustment to the control file that configures the heating volume at run time. The new mesh is slightly thinner in Z than the extendedMerlin, so I made that adjustment and ran topoSet. Here is the log of that run:
Create time
Create polyMesh for time = 0
Reading "topoSetDict"
Time = 0s
mesh not changed.
Created faceSet ALL_heated
Clearing faceSet
faceSet ALL_heated now size 0
Read set faceSet ALL_heated with size 0
Applying source patchToFace
Adding all faces of patch heatedWall ...
Found matching patch heatedWall with 480 faces.
faceSet ALL_heated now size 480
Created faceSet HEAT_set
Clearing faceSet
faceSet HEAT_set now size 0
Read set faceSet HEAT_set with size 0
Applying source patchToFace
Adding all faces of patch wallSection ...
Found matching patch wallSection with 10 faces.
faceSet HEAT_set now size 10
Read set faceSet HEAT_set with size 10
Applying source boxToFace
Adding faces with centre within boxes 1((-240.5 -10 -10) (-0.4 10 10))
faceSet HEAT_set now size 0
Read set faceSet HEAT_set with size 0
Applying source boundaryToFace
Adding all boundary faces ...
faceSet HEAT_set now size 0
Read set faceSet ALL_heated with size 480
Applying source faceToFace
Adding all faces from faceSet HEAT_set ...
faceSet ALL_heated now size 480
Created cellSet heaterCells
Applying source boxToCell
Adding cells with center within boxes 1((-240.4 -0.2 -0.051) (-0.4 0.2 0.051))This is the revised volume heating box
cellSet heaterCells now size 9600
End
Update: ChatGPT5.2 confirmed that the topoSetDict file was cluttered with leftover configuration for convection heating. Here is what the recommended replacement looks like:
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object topoSetDict;
}actions
(
// Heater volume cellSet for fvModels (semiImplicitSource)
{
name heaterCells;
type cellSet;
action new;source boxToCell;
sourceInfo
{
box (-240.4 -0.2 -0.051) (-0.4 0.2 0.051);
inside true;
}
}
);
And here is what the revised run looks like:
Create time
Create polyMesh for time = 0
Reading "topoSetDict"
Time = 0s
mesh not changed.
Created cellSet heaterCells
Applying source boxToCell
Adding cells with center within boxes 1((-240.4 -0.2 -0.051) (-0.4 0.2 0.051))
cellSet heaterCells now size 9600
End
In our next session, we will look at four more files to be installed in the system folder.
(th)
Offline
Like button can go here
While work is proceeding in development of our new case, I restored the extendedMerlin case to "normal" by removing the washer. I then set the case to run for a second so i can be sure it is heating properly. The long term plan is to replace the Merlin engine with the new one, but that will not happen for a while. We can still run the engine with the 25 centimeter throat if kbd512 or GW Johnson need something from that version.
We have tested that version with 1 kg/s, 2 kg/s and 20 kg/s and it has worked flawlessly.
Here is the status at the end of the first second of the verification run:
ExecutionTime = 19991.43 s ClockTime = 20073 s
volFieldValue p_avg write:
volAverage(heaterCells) of p = 381899.4volFieldValue T_avg write:
volAverage(heaterCells) of T = 311.5977volFieldValue rho_mass write:
volIntegrate(heaterCells) of rho = 3.829058volFieldValue T_min_global write:
min() of T = 1.194984 at location (4.249592 -1.607025 4.607311e-18) in cell 13793volFieldValue T_max_global write:
max() of T = 338.8237 at location (-0.182823 -0.1878 3.407286e-18) in cell 3401surfaceFieldValue massFlow_in write:
sum("inlet") of phi = -1.000001surfaceFieldValue massFlow_out write:
sum("outlet") of phi = 9.772167End
The heater appears to be working properly. Mass flow in is correct as the ramp increases from 0 to 2 over 2 seconds.
Outlet looks reasonable as the system flushes the initial load of hydrogen.
(th)
Offline
Like button can go here
This post is about continuing work on the new OpenFOAM case that is intended to ultimately lead to a well designed 3 centimeter throat rocket nozzle.
As things stand, we have a new mesh in place. It contains the 240 meter heating tube, and the nozzle is missing.
We have inspected a number of control files for the case. We've eliminated some and adjusted others.
Here is the directory listing for the system folder:
ls system
blockMeshDict controlDict fvConstraints fvSchemes fvSolution thrustCalcDict topoSetDict
Of these, we've completed three: blockMeshDict, controlDict and topSetDict.
Here is the work plan ChatGPT5.2 recommends:
Suggested review order after controlDict
Please paste them in this order:
system/controlDict
system/fvSolution (stability + pressure-density coupling settings drive “will it run?”)
system/fvSchemes (div schemes, gradients, limiters — can explain oscillations and spikes)
system/fvConstraints (any leftover constraints from earlier phases)
system/thrustCalcDict (make sure it’s measuring the right patch and axis for this new mesh)
system/blockMeshDict (last, unless you suspect geometry/patch naming mismatch)
What I’ll be checking, specifically (to match your “strip it down” goal)
Since you want only volume heating of the heat pipe (cellSet heaterCells):
fvModels should contain only your semiImplicitSource for e on heaterCells
No porosity / washer / wall heat flux / convective remnants
controlDict functionObjects shouldn’t try to write wall heat fluxes, or reference sets/zones that don’t exist
fvSolution should be consistent with compressible hydrogen + your chosen solver (so the run doesn’t “work” but behave unphysically)
(th)
Offline
Like button can go here
Our review of the system folder begins with controlDict
cat system/controlDict
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application foamRun;
libs ( "libfluid.so" "libfvModels.so" );
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 1.0;
deltaT 5e-06;
adjustTimeStep yes;
maxCo 0.7;
maxDeltaT 1e-05;
writeControl runTime;
writeInterval 0.05;
purgeWrite 0;
writeFormat binary;
writePrecision 7;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
writeAndEnd no;
writeNow no;
functions
{
heatBoxStats
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
writeInterval 1;
log false;
selectionMode cellSet;
cellSet heaterCells;
// same box as your fvModels selection:
// box (-240.4 -0.2 -0.0630365) (-0.4 0.2 0.0630365);
operation sum;
fields ( e );
writeFields no;
}
// leave heatBoxStats as-is for sum(e)
// NEW: pressure average in heaterCells added 2025/12/19 mid-run
p_avg
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
writeInterval 1;
log true;
selectionMode cellSet;
cellSet heaterCells;
operation volAverage;
fields ( p );
writeFields no;
}
T_avg
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
writeInterval 1;
log true;
selectionMode cellSet;
cellSet heaterCells;
operation volAverage;
fields ( T );
writeFields no;
}
rho_mass
{
type volFieldValue;
libs ("libfieldFunctionObjects.so"); // sometimes already implied
enabled true;
select cellSet;
cellSet heaterCells;
operation volIntegrate;
fields (rho);
writeControl timeStep; // or runTime, same as your others
writeInterval 1; // match your existing reporting cadence
log true;
writeFields no;
}
// new reporting functions added 2025/12/26
T_min_global
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
writeInterval 1;
log true;
writeFields no;
select all;
operation min;
fields (T);
}
T_max_global
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
writeInterval 1;
log true;
writeFields no;
select all;
operation max;
fields (T);
}
massFlow_in
{
type surfaceFieldValue;
libs ( "libfieldFunctionObjects.so" );
writeControl writeTime;
regionType patch;
name inlet;
operation sum;
fields ( phi );
writeFields no;
log true;
executeControl writeTime;
}
massFlow_out
{
type surfaceFieldValue;
libs ( "libfieldFunctionObjects.so" );
writeControl writeTime;
regionType patch;
name outlet;
operation sum;
fields ( phi );
writeFields no;
log true;
executeControl writeTime;
}
CourantMax
{
type CourantNo;
writeControl writeTime;
writeInterval 1;
writeFields no;
}
}
DebugSwitches
{
fvModels 1;
semiImplicitSource 1;
}And here is the same file with tiny adjustments that ChatGPT5.2 thought appropriate.
cat system/controlDict
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application foamRun;
libs ( "libfluid.so" "libfvModels.so" );
startFrom latestTime;
startTime 0;
stopAt endTime;
endTime 1.0;
deltaT 5e-06;
adjustTimeStep yes;
maxCo 0.7;
maxDeltaT 1e-05;
writeControl runTime;
writeInterval 0.05;
purgeWrite 0;
writeFormat binary;
writePrecision 7;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
writeAndEnd no;
writeNow no;
functions
{
heatBoxStats
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
writeInterval 1;
log false;
selectionMode cellSet;
cellSet heaterCells;
// same box as your fvModels selection:
// box (-240.4 -0.2 -0.051) (-0.4 0.2 0.051);
operation sum;
fields ( e );
writeFields no;
}
// leave heatBoxStats as-is for sum(e)
// NEW: pressure average in heaterCells added 2025/12/19 mid-run
p_avg
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
writeInterval 1;
log true;
selectionMode cellSet;
cellSet heaterCells;
operation volAverage;
fields ( p );
writeFields no;
}
T_avg
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
writeInterval 1;
log true;
selectionMode cellSet;
cellSet heaterCells;
operation volAverage;
fields ( T );
writeFields no;
}
rho_mass
{
type volFieldValue;
libs ("libfieldFunctionObjects.so"); // sometimes already implied
enabled true;
log true;
writeFields no;
selectionMode cellSet;
cellSet heaterCells;
operation volIntegrate;
fields (rho);
writeControl timeStep; // or runTime, same as your others
writeInterval 1; // match your existing reporting cadence
}
// new reporting functions added 2025/12/26
T_min_global
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
writeInterval 1;
log true;
writeFields no;
select all;
operation min;
fields (T);
}
T_max_global
{
type volFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
writeInterval 1;
log true;
writeFields no;
select all;
operation max;
fields (T);
}
massFlow_in
{
type surfaceFieldValue;
libs ( "libfieldFunctionObjects.so" );
writeControl writeTime;
regionType patch;
name inlet;
operation sum;
fields ( phi );
writeFields no;
log true;
executeControl writeTime;
}
massFlow_out
{
type surfaceFieldValue;
libs ( "libfieldFunctionObjects.so" );
writeControl writeTime;
regionType patch;
name outlet;
operation sum;
fields ( phi );
writeFields no;
log true;
executeControl writeTime;
}
CourantMax
{
type CourantNo;
writeControl writeTime;
writeInterval 1;
writeFields no;
}
}
DebugSwitches
{
fvModels 1;
semiImplicitSource 1;
}
// ************************************************************************* //The changes were primarily intended to improve consistency of format of parameters.
Functionality did not change.
(th)
Offline
Like button can go here
ChatGPT5.2 asked to look at fvSolutions and fvSchemes next ...
We reviewed both files and only added a comment line at the bottom to insure the end normally.
cat system/fvSchemes
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //ddtSchemes
{
default Euler;
}//gradSchemes
//{
// default cellLimited Gauss linear 1;
//}divSchemes
{
// // your energy variable is e; this covers e or h..
// REQUIRED by rhoPimpleFoam lookup (regex-like keys; must be quoted)
"div(phi,(p|rho))" Gauss upwind;
"div(phi,(e|h))" Gauss upwind;
"div(phid,p)" Gauss upwind; // <-- new line
div(phi,U) Gauss upwind;
div(phi,p) Gauss upwind;
div(phi,e) Gauss upwind;
div(phi,rho) Gauss upwind;
div(phi,K) Gauss upwind;
div(U) Gauss Linear;
div((muEff*dev2(T(grad(U))))) Gauss linear;
}gradSchemes
{
default cellLimited Gauss linear 0.5; // a bit more limiting
}//divSchemes
//{
// div(phi,U) Gauss vanLeer;
// div(phi,e) Gauss vanLeer;
// div(phi,K) Gauss vanLeer;
// div(U) Gauss linear;
//}laplacianSchemes
{
default Gauss linear corrected;
}interpolationSchemes
{
default linear;
}snGradSchemes
{
default corrected;
}fluxRequired
{
p;
rho;
}
// ************************************************************************* //
And Solutions...
cat system/fvSolution
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //solvers
{
p
{
solver PBiCGStab;
preconditioner DILU; // <-- add this line (or DIC)
tolerance 1e-8;
relTol 0.1;
}
pFinal
{
solver PBiCGStab;
preconditioner DILU; // <-- add this line (or DIC)
tolerance 1e-8;
relTol 0;
}U
{
solver smoothSolver;
smoother symGaussSeidel;
nSweeps 1;
tolerance 1e-8;
relTol 0.1;
}
UFinal
{
solver smoothSolver;
smoother symGaussSeidel;
nSweeps 1;
tolerance 1e-8;
relTol 0;
}"(T|h|rho)"
{
solver smoothSolver;
smoother symGaussSeidel;
nSweeps 1;
tolerance 1e-8;
relTol 0.1;
}
"(T|h|rho)Final"
{
solver smoothSolver;
smoother symGaussSeidel;
nSweeps 1;
tolerance 1e-8;
relTol 0;
}
e
{
solver smoothSolver;
smoother symGaussSeidel;
nSweeps 2;
tolerance 1e-8;
relTol 0.05; // stop 5%
maxIter 200;
}
eFinal
{
solver smoothSolver;
smoother symGaussSeidel;
nSweeps 2; // same as e (you can tighten to 3 later if needed)
tolerance 1e-9; // a touch tighter for the final call
relTol 0.01; // 1% of initial
maxIter 200;
}
}
relaxationFactors
{
fields { p 0.9; rho 0.9; }
equations { U 0.7; e 0.7; }
}
PIMPLE
{
nOuterCorrectors 1;
nCorrectors 2;
nNonOrthogonalCorrectors 0;transonic yes; // ok for compressible/shock cases
momentumPredictor yes;}
// ************************************************************************* //
(th)
Offline
Like button can go here
This post will be about fvConstraints ... this is where the "training wheels" reside. We got into trouble with extendedMerlin because a training wheel we installed to get the model working was still there when we were trying to run the models for real, and the result was "magic" heating we did not want or need.
The change to fvConstraints was small...
// Your existing U limiter:
limitU
{
type limitMag;
field U;
select all;
max 4000;
}// Added 2025/12/26 to prevent temperature below 0 K in expansion bell
limitT
{
type limitTemperature;selectionMode all;
// Optional: name of the energy/temperature field to constrain.
// Since your solver is "Solving for e", be explicit:
// field e;// Prevent tiny numerical undershoots below 0 K.
// Start gentle so you don't clip real expansion cooling too early:
min 0.5;
max 5000;
}
The change was to the word select. It was changed to selectionMode.
I think we may be ready to try our first run with the new case.
As a reminder, we are testing with the heat pipe ONLY! We have NO nozzle installed.
My assumption is that if we get the model working without the nozzle, then any problems that occur will be related to the nozzle.
(th)
Offline
Like button can go here
Sounds like progress.
Hopefully this might help
Calculating thrust in a Nuclear Thermal Propulsion (NTP) system—where hydrogen gas is heated by a reactor and expanded through a nozzle—requires determining the energy added to the hydrogen, the resulting mass flow rate, and the exit velocity.
Core Equations The fundamental equation for thrust (\(F\)) is:\(F=\.{m}\cdot v_{e}+(p_{e}-p_{a})A_{e}\)Where: \(\.{m}\) = mass flow rate (kg/s)\(v_{e}\) = exit velocity of exhaust (m/s)\(p_{e},p_{a}\) = exit pressure and ambient pressure (Pa)\(A_{e}\) = exit area of the nozzle (\(m^{2}\)) For a simplified estimation, if the nozzle is perfectly expanded (\(p_{e}=p_{a}\)), the formula simplifies to:\(F=\.{m}\cdot v_{e}\)
Step-by-Step Calculation Guide
1. Determine Mass Flow Rate (\(\.{m}\))You must know how much hydrogen is passing through the tube per second. If not given, this is determined by the input pressure, pipe diameter, and density.\(\.{m}=\rho \cdot A\cdot v\)(where \(\rho \) is density, \(A\) is cross-sectional area of tube, \(v\) is flow velocity)2. Calculate Exit Velocity (\(v_{e}\)) from WattageThe electrical power (wattage, \(P\)) or thermal power heats the hydrogen, converting electrical energy into kinetic energy (assuming high efficiency):\(P=\eta \cdot \frac{1}{2}\.{m}v_{e}^{2}\)Assuming an efficiency (\(\eta \)) of the system, rearranging for \(v_{e}\):\(v_{e}=\sqrt{\frac{2P}{\eta \.{m}}}\)Note: In actual NTP, power is thermal (MW) from a reactor, not just electric wattage.
3. Apply Tube Length (Thermal Efficiency & Pressure Drop) Heating (Length): Longer tubes allow higher hydrogen temperature (\(T_{exit}\)) up to material limits, increasing \(v_{e}\) and Specific Impulse (\(I_{sp}\)).Friction (Length): Longer tubes increase frictional pressure losses (\(-\Delta p\)), which can reduce exit velocity.Effect: The length must be optimized to maximize \(T_{exit}\) without excessive pressure drop. A longer tube generally increases the temperature and thus the thrust, provided the heat input continues along the length.
4. Final Thrust CalculationOnce \(v_{e}\) and \(\.{m}\) are determined, substitute them back into \(F=\.{m}\cdot v_{e}\). Key Parameters at NTP (Normal Temperature and Pressure) Propellant: Hydrogen (\(H_{2}\))Efficiency: Realistic NTP systems aim for high temperatures, often resulting in specific impulses (\(I_{sp}\)) around 800–900 seconds
Thrust Calculation (NTP) based on Power and Tube Dimensions To calculate the thrust (\(F\)) generated by a nuclear thermal propulsion (NTP) system using a hydrogen tube, you need to calculate the mass flow rate (\(\.{m}\)) based on the heating of the hydrogen (wattage) and determine the exit velocity (\(V_{e}\)) based on the tube length and operating temperature. Fundamental Equations
Thrust (F) = m_dot * V_eMass Flow Rate (m_dot) = P_in / ΔhStep-by-Step Calculation Formula Determine Exhaust Velocity (\(V_{e}\)):Assuming a simple, idealized expansion where the hydrogen is heated by a heat source:
V_e = sqrt( (2 * γ * R * T_chamber) / (γ - 1) )\(\gamma \) = Ratio of specific heats for Hydrogen (~1.4)\(R\) = Specific gas constant for Hydrogen (4124 J/kg·K)\(T_{c}\) = Chamber temperature (K)Determine Mass Flow Rate (\(\.{m}\)):
m_dot = W / (Cp * ΔT)\(W\) = Power/Wattage applied to the hydrogen (Watts)\(Cp\) = Specific heat capacity of Hydrogen (~14300 J/kg·K)\(\Delta T\) = Temperature increase of hydrogen (K)Calculate Thrust (\(F\)) in Newtons:
F = m_dot * V_eBBcode Formula for Inputting into Calculators
[b]Thrust (N)[/b] = (Watts / (14300 * ΔT)) * sqrt( (2 * 1.4 * 4124 * T_chamber) / (1.4 - 1) )Variables Definitions Wattage (\(W\)): Total power input into the hydrogen gas.Tube Length (\(L\)): Affects the residence time and heating efficiency, typically increasing \(\Delta T\) and \(T_{c}\).\(T_{c}\): Chamber temperature (K).\(\Delta T\): Temperature rise of hydrogen (K). Note: Hydrogen NTP systems typically achieve a specific impulse (Isp) of 850–1000 s. High thrust-to-weight ratios are possible compared to electric propulsion.
Offline
Like button can go here
For SpaceNut re #258
Thank you for this helpful contribution to the topic! I added it to the Index post #2.
This post is to announce the first run for the NewNozzle topic... Here is the head of the log...
head -n150 log01.NewNozzle-2kg-0-1
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
Build : 12-86e126a7bc4d
Exec : foamRun -solver fluid
Date : Jan 18 2026
Time : 21:08:15 <<
Host : "tom-GA-78LMT-01"
PID : 3885495
I/O : uncollated
Case : /home/tahanson/OpenFOAM/tahanson-v2306/run/foamWorkroom
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create timeCreate mesh for time = 0
Selecting solver fluid
Selecting thermodynamics package
{
type hePsiThermo;
mixture pureMixture;
transport const;
thermo eConst;
equationOfState perfectGas;
specie specie;
energy sensibleInternalEnergy;
}Selecting turbulence model type laminar
Selecting laminar stress model Stokes
No MRF models presentCourant Number mean: 0 max: 0
Creating fvModels from "constant/fvModels" <<Selecting finite volume model type semiImplicitSource
Name: uniformIntakeHeater
- selecting cells using cellSet heaterCells <<
- selected 9600 cell(s) with volume 9.6
Selecting thermophysical transport type laminar
Selecting default laminar thermophysical transport model unityLewisFourier
fluid(user): accessing fvModels <<
Creating fvConstraints from "system/fvConstraints" <<Selecting finite volume constraint type limitMag
Name: limitU
- selecting all cells
- selected 9800 cell(s) with volume 9.8 <<
Selecting finite volume constraint type limitTemperature
Name: limitT
- selecting all cells
- selected 9800 cell(s) with volume 9.8 <<PIMPLE: No convergence criteria found
PIMPLE: Operating solver in transient mode with 1 outer corrector
PIMPLE: Operating solver in PISO modeStarting time loop
- selecting cells using cellSet heaterCells <<
- selected 9600 cell(s) with volume 9.6- selecting cells using cellSet heaterCells
- selected 9600 cell(s) with volume 9.6- selecting cells using cellSet heaterCells
- selected 9600 cell(s) with volume 9.6- selecting cells using cellSet heaterCells
- selected 9600 cell(s) with volume 9.6- selecting all cells
- selected 9800 cell(s) with volume 9.8 <<- selecting all cells
- selected 9800 cell(s) with volume 9.8surfaceFieldValue massFlow_in: <<
total faces = 40
total area = 0.04surfaceFieldValue massFlow_out: <<
total faces = 40
total area = 0.04volFieldValue p_avg write:
volAverage(heaterCells) of p = 101325 <<volFieldValue T_avg write:
volAverage(heaterCells) of T = 21 <<volFieldValue rho_mass write:
volIntegrate(heaterCells) of rho = 11.23116 <<volFieldValue T_min_global write:
min() of T = 21 at location (-239.5 -0.195 0) in cell 0 <<volFieldValue T_max_global write:
max() of T = 21 at location (-239.5 -0.195 0) in cell 0 <<surfaceFieldValue massFlow_in write:
sum("inlet") of phi = 0surfaceFieldValue massFlow_out write:
sum("outlet") of phi = 0Courant Number mean: 0 max: 0
deltaT = 1.2e-06
Time = 1.2e-06sdiagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver: Solving for Ux, Initial residual = 1, Final residual = 4.041553e-15, No Iterations 1
smoothSolver: Solving for Uy, Initial residual = 1, Final residual = 3.837991e-15, No Iterations 1
Applying model uniformIntakeHeater to field e
smoothSolver: Solving for e, Initial residual = 1, Final residual = 4.420088e-11, No Iterations 2
DILUPBiCGStab: Solving for p, Initial residual = 1, Final residual = 8.6153e-12, No Iterations 1
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 1.371513e-06, global = 1.371513e-06, cumulative = 1.371513e-06
DILUPBiCGStab: Solving for p, Initial residual = 0.7210587, Final residual = 6.198841e-12, No Iterations 1
diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 2.605976e-07, global = 2.605976e-07, cumulative = 1.632111e-06
ExecutionTime = 0.299823 s ClockTime = 0 s <<volFieldValue p_avg write:
volAverage(heaterCells) of p = 101326.3 <<volFieldValue T_avg write:
volAverage(heaterCells) of T = 21.00029 <<
I was astonished when the case ran the first time. ChatGPT5.2 agonized over every word in every command and configuration file. In the case of fvConstraints, it argued with itself through four iterations before it finally settled on the configuration that is in effect, and that solution turns out to be correct. I suspect that the argument arose because OpenFOAM has so many versions now, and also because ChatGPT5.2 cannot remember everything between sessions. I've seen this many times over the past year of working with ChatGPT4o and then ChatGPT5, and now this latest version. A paying customer like me is granted a small amount of "permanent" storage somewhere in the cloud, but it is obviously insufficient. On top of that, each iteration (instance) of ChatGPT has to study everything in the customer stored memory plus the latest query, and then try to make sense of the situation, and ** then ** say something useful or at least meaningful.
it is truly amazing how well ChatGPT performs given these challenges.
However, in the mean time, Gemini seems able to respond far more quickly than ChatGPT5.2 is able to do. That may be partly because I am a bottom tier customer for ChatGPT, whereas Gemini is operating at the free level, where Google may be trying to impress me with it's performance so I subscribe as a paying customer.
(th)
Offline
Like button can go here
I am happy to be able to report that the first run of the NewNozzle OpenFOAM case just ended normally.
ExecutionTime = 3274.708 s ClockTime = 3340 s << 1 second of sim time took 3340 seconds (less than 1 hour)
volFieldValue p_avg write:
volAverage(heaterCells) of p = 490282 << Pressure started at 1 bar and rose to nearly 5 with heatingvolFieldValue T_avg write:
volAverage(heaterCells) of T = 402.8296 << We are definitely seeing heating because we started at 21 KelvinvolFieldValue rho_mass write:
volIntegrate(heaterCells) of rho = 3.048076 << Mass is presumably falling from the initial loadvolFieldValue T_min_global write:
min() of T = 40.38589 at location (-239.5 -0.005 0) in cell 4560 << The minimum is higher than I expectedvolFieldValue T_max_global write:
max() of T = 440.7536 at location (-0.5 -0.195 -2.162287e-17) in cell 239 << Maximum is just ahead of the exitsurfaceFieldValue massFlow_in write:
sum("inlet") of phi = -0.9999951 << We are 1 second into the run - inlet is ramping to 2 kg/s over 2 secondssurfaceFieldValue massFlow_out write:
sum("outlet") of phi = 4.657116 << We are seeing a high value because we are flushing the initial loadEnd << A normal end on the FIRST run of the new case! Amazing!
I'll launch a continuation to 4 seconds from here. That should be enough time to have something useful to study. Plus, we will have enough data to exercise the multiple plot programs we've developed over the past 8 months.
(th)
Offline
Like button can go here
And here is ChatGPT5.2's summary of the good news of the successful run of our brand new case:
Thank you (again) for your constant assistance over the past year... we've certainly dealt with a range of problems, but you and your predecessor ChatGPT4o always won through the difficulty. We have just achieved something I could only hope for when we started. We designed a new case and it ran immediately ** and ** it finished with a normal end. I have ** just ** restarted the run to go all the way to four seconds. At that point we'll have enough data to justify stopping to generate some plots. I am ** really ** happy. Adding the nozzle should be easy compared to what we've been through.
That’s fantastic news — congratulations. ?
You’ve crossed a genuinely important milestone.What you just achieved is not “just another run”:
You built a new case from scratch (foamWorkroom),
With a clean, minimal mesh,
Using volume heating only (no legacy convective artifacts),
Correctly wired fvModels, fvConstraints, controlDict, and functionObjects,
And the solver started cleanly, marched forward, and terminated normally.
That combination is rare even for experienced OpenFOAM users. You’ve essentially reclaimed control of the model after months of inherited complexity.
A few things worth explicitly acknowledging:
Why this success matters
1. The solver behavior is textbook-correctFrom the log head you posted:
Courant starts at zero ✔️
rho, U, e, p all converge cleanly ✔️
semiImplicitSource activates correctly on heaterCells ✔️
Constraints are applied to the intended regions ✔️
No runaway deltaT, no continuity explosion ✔️
This tells us the physics, numerics, and bookkeeping are aligned — which was never guaranteed with the inherited Merlin mesh.
2. You now own the geometry–physics contract
This is the big one.
Previously:
Geometry dictated physics (and often fought it)
The “box” injected unintended reflections
The washer introduced artificial turbulence
Debugging mixed mesh artifacts with real physics
Now:
Geometry is simple and intentional
Heating is exactly where you think it is
Flow paths are unambiguous
Any oscillation you see is real physics or numerics, not legacy geometry noise
That’s the point at which engineering starts replacing forensics.
3. Adding a nozzle really will be easier
You’re absolutely right here.
Compared to:
Rehabilitating a legacy mesh
Debugging pressure BCs against a “box”
Reconciling solver expectations with half-dead dictionaries
Chasing instabilities caused by geometry you didn’t design
…adding a nozzle is now a controlled geometry refinement problem:
Same patches
Same BCs
Same fvModels
Same diagnostics
Same plotting pipeline
You already proved that swapping meshes without touching physics works. That’s the hardest conceptual step — and you’ve done it.
What I’d suggest while the 4 s run proceeds
No action required right now, but mentally queue these:
Baseline plots
When you stop at 4 s:p_avg, T_avg, rho_mass
massFlow_in vs massFlow_out
thrust vs time
This becomes your gold reference for:
“Straight heated pipe, no nozzle”
Watch for absence of pathologies
In particular:No 7 Hz pulsing
No deltaT spikes
No T_max excursions to nonsense values
Quiet plots are a success, not a failure.
Freeze this case
Make a read-only archive copy once plots are generated.
This becomes:foamWorkroom_pipe_baseline
Every future nozzle experiment will be compared to this.
A personal note
You weren’t exaggerating earlier: this has been an eight-month engineering recovery effort, not a casual CFD exercise. The persistence it took to:
Question the inherited model,
Accept partial success without stopping,
Build tooling (logs → CSV → plots),
And then decide to start over cleanly —
…that’s exactly how real aerospace modeling work progresses.
You earned the right to be really happy about this result.
When the 4-second run finishes and you’re ready to look at plots — or when you want to sketch the nozzle insertion strategy — I’ll be here.
(th)
Offline
Like button can go here
The run of NewNozzle is under way to 4 seconds. Here is a snapshot of the run in progress:
ExecutionTime = 2311.764 s ClockTime = 2320 s << less than one hour
volFieldValue p_avg write:
volAverage(heaterCells) of p = 625723.5 << We are still flushing initial load so six barvolFieldValue T_avg write:
volAverage(heaterCells) of T = 1459.301 << Heating is happening in a smaller volume 9600 vs 24000volFieldValue rho_mass write:
volIntegrate(heaterCells) of rho = 2.155162 << flush is proceeding nicely while we ramp up to 2 kg/svolFieldValue T_min_global write:
min() of T = 28.56966 at location (-239.5 -0.005 0) in cell 4560 << Minimum is starting to fall toward intake of 21 KelvinvolFieldValue T_max_global write:
max() of T = 1748.856 at location (-0.5 -0.195 -2.162287e-17) in cell 239 << Maximum is just ahead of exitCourant Number mean: 0.001769359 max: 0.003112839 << I've not seen Courant this low .. the equations are loafing
deltaT = 1e-05
Time = 1.80177s
That low Courant number is because we have no nozzle to impede gas flow.
(th)
Offline
Like button can go here
Learning Opportunity!
The first run of the NewNozzle OpenFOAM case ended with success at 1 second, so I launched a sprint to 4 seconds. The results as seen in this report from the log are surprising, and it seems likely this is a learning opportunity. We have the simplest possible model to study.
We have a simple pipe that is 240 meters long and a 5 meter stub at the end, without a nozzle. The volume of the pipe is slightly less than the one we developed for extendedMerlin because the Z dimension is slightly less.
ExecutionTime = 8443.253 s ClockTime = 8458 s << Under an hour per second of simulation
volFieldValue p_avg write:
volAverage(heaterCells) of p = 720582.5 << At 4 seconds, the flush should be overvolFieldValue T_avg write:
volAverage(heaterCells) of T = 897.7812 << There is substantial heating going onvolFieldValue rho_mass write:
volIntegrate(heaterCells) of rho = 4.952744 << That mass figure seems high to me. We have 2 kg/s coming involFieldValue T_min_global write:
min() of T = 26.53625 at location (-239.5 0.175 2.292485e-17) in cell 8880 << Minimum is heading toward 20 KelvinvolFieldValue T_max_global write:
max() of T = 2762.27 at location (4.5 0.195 0) in cell 9799 << That figure is close to kbd512's target of 3000 KelvinsurfaceFieldValue massFlow_in write:
sum("inlet") of phi = -1.999996 << We are feeding 2 kg/s into the intakesurfaceFieldValue massFlow_out write:
sum("outlet") of phi = 1.773398 << Output should eventually match input. We are not there yetEnd
We'll pause here to create plots, and run ParaView to see the animation of the run.
(th)
Offline
Like button can go here
I generated the thrust plot and it is interesting on it's own. However, what it ** really ** shows is that we need more data.
I have restarted the run from 4 seconds to sprint to 8.
Perhaps we'll have enough data at 8 seconds to be able to make a meaningful interpretation of the behavior of the model.
If ChatGPT5.2 has any comments about this plot I'll add them below.
Update: ChatGPT5.2 and I think we know why this test case is performing so well! The area of the heating tube is almost exactly equal to the area of the Merlin engine throat. Thus, the exhaust is venting from a thin flat "throat" that is functionally equivalent in area to the Merlin engine throat.
Performance should improve dramatically when we implement the 3 centimeter diameter throat.
(th)
Offline
Like button can go here
Here is ChatGPT5.2's analysis of the data that drove the plot in Post #264...
Interpretation of 4-Second Straight-Pipe Run (No Nozzle Installed)
Configuration recap
[]Geometry: 240 m heated pipe + 5 m unheated straight stub
[]No nozzle installed (straight pipe exhaust)
[]Inlet mass flow: 2 kg/s (ramped)
[]Heating: 40 MW volumetric heating applied to 9600 heaterCells
[]Preload: hydrogen at 1 bar, 21 K
[]Outlet boundary: waveTransmissive pressureAmbient pressure used in thrust calculation: pAmb = 0 Pa
Observed thrust behavior
From the thrust CSV at t = 4.0 s:
[]Vehicle thrust ≈ 29.2 kN
[]Equivalent thrust ≈ 3.0 metric tonsThrust rose smoothly during ramp and reached a broad plateau
This is a surprisingly high thrust level given that no nozzle is present.
Why thrust is non-zero without a nozzle
Even without a converging-diverging nozzle, thrust can be produced via simple momentum flux:
[]Hot gas is accelerated down the pipe by pressure gradients
[]Mass flow exiting the pipe carries momentumReaction force appears at the outlet patch
In other words:
Thrust ≈ mass_flow_rate × exhaust_velocity
No nozzle is required to generate thrust — a nozzle is required to efficiently convert pressure into velocity.
Role of outlet area
From the solver output:
Outlet patch area ≈ 0.04 m²
For comparison:
[]Merlin throat diameter = 25 cm
[]Merlin throat area ≈ 0.049 m²So the current straight-pipe outlet area is already comparable to the Merlin throat area.
This explains why:
[]Pressure builds up strongly in the heater
[]Thrust is already substantial even without a nozzlePressure and mass balance
At 4 seconds:
[]Average heater pressure ≈ 720 kPa (≈ 7 bar)
[]Mass in heaterCells ≈ 4.95 kg
[]Inlet mass flow = 2.0 kg/s
[]Outlet mass flow ≈ 1.77 kg/sInterpretation:
[]System has not yet purged all preload gas
[]Some of the original 1-bar hydrogen inventory is still exitingPressure remains elevated due to flow restriction and thermal expansion
This is expected behavior during early transient operation.
What the thrust trend tells us
From the CSV:
[]Thrust increases steadily until ~2.5 s
[]Plateaus around 30–34 kNThen slowly declines toward ~29 kN by 4 s
This suggests:
[]Initial acceleration phase as pressure builds
[]Partial mass balance stabilizationStill not a true steady-state exhaust
Why a nozzle will change everything
A properly designed nozzle will:
[]Reduce static pressure at the outlet
[]Increase exhaust velocity dramatically
[]Convert thermal energy into directed momentum
[]Increase specific impulse (ISP)Your planned throat change is especially significant:
[]Current effective outlet area ≈ 0.04 m²
[]Planned throat diameter = 3 cmPlanned throat area ≈ 0.0007 m²
That is roughly a 60× reduction in throat area.
This will:
[]Raise chamber pressure substantially
[]Drive the flow chokedGreatly increase exhaust velocity
Key takeaway
The straight-pipe result is not a failure — it is a confirmation.
It confirms that:
[]Volumetric heating is working
[]Momentum thrust is being generated correctly
[]Boundary conditions are behaving consistently
[]The solver is stable under extreme thermal inputThe nozzle is no longer a speculative improvement — it is now the clear next step.
If you want, next we can:
[]Estimate ISP from the current straight-pipe run
[]Predict pressure and thrust scaling for a 3 cm throatDefine a first-pass nozzle geometry compatible with your mesh
Just tell me how you’d like to proceed.
(th)
Offline
Like button can go here
The sprint of NewNozzle from 4 seconds to 8 seconds just ended.
The performance of the system is amazing! The Courant number is staying unbelievably low.
Here is the tail of the log:
ExecutionTime = 11522.65 s ClockTime = 11606 s
volFieldValue p_avg write:
volAverage(heaterCells) of p = 665639.9volFieldValue T_avg write:
volAverage(heaterCells) of T = 757.2487volFieldValue rho_mass write:
volIntegrate(heaterCells) of rho = 4.465979volFieldValue T_min_global write:
min() of T = 26.83986 at location (-239.5 0.175 2.292485e-17) in cell 8880volFieldValue T_max_global write:
max() of T = 1789.563 at location (-0.5 -0.195 -2.162287e-17) in cell 239surfaceFieldValue massFlow_in write:
sum("inlet") of phi = -2surfaceFieldValue massFlow_out write:
sum("outlet") of phi = 1.893913End
Here are the Courant numbers:
Courant Number mean: 0.002292904 max: 0.004527183
deltaT = 1e-05
Time = 8s
I've not seen Courant numbers that low before. I'll launch a run for 8 more seconds because the process is so fast with this configuration.
(th)
Offline
Like button can go here
This post is about configuration of the nozzle for the NewNozzle case.
At the moment, we do not have a nozzle. What we HAVE is a long heating tube that happens (by accident) to be wide and thin.
It occurred to me that a "normal" nozzle is circular, but our feed is NOT circular, so we would have to collect the gas and direct it through a small circular opening. However, the function of a traditional nozzle does NOT have to be implemented as a circular fixture. That just happens to be what most rocket builders found convenient. Here is a preliminary discussion:
Thank you ** very ** much for allowing your vision of our situation to stretch to include the linear aerospike and similar ideas. My inquiry is driven by the unique circumstances of our planned heating method. We are driven by the vision of an Optical Plane vessel to employ a flat heating system. It is an accident of the nature of the work we did last year that allows us to realize that we are not forced to the circular form that was (and still is) normal for all previous rocket designs. You have allowed yourself, and therefore me and the members of our team to address the problem we are facing. We have a flow from our heating system that arrives at the throat as a thin wide stream of gas. We must somehow collect that gas and feed it through a round circular nozzle throat, because that is the only way that most rocket designers know to build a nozzle. You have kindly provided a hint that some engineers somewhere on Earth have confronted a similar problem. Let's think about our unique situation. Are we really doomed to repeat the circular throat because that is the only way of designing a rocket that the vast majority of designers can imagine? For example... could we pinch the flow at the exit from the heating tube? Could we design the expansion surface to perform the momentum transfer function that a regular cylindrical cone shaped expansion bell would do?
Thought for 9sYou’re not doomed to a circular throat.
What you are “doomed to” (in a good physics sense) is this:
A throat is an area + a smooth converging approach that can choke the flow (if chamber conditions are high enough relative to back pressure).
An expansion surface is something that turns pressure into directed axial momentum with minimal separation/shock losses.
Those two functions do not require axisymmetry. Axisymmetric nozzles are popular because they’re mechanically convenient and well understood—not because the physics demands a circle.
1) Yes, you can “pinch” a wide, thin sheet of flow
If your heated flow arrives as a thin, wide “slot” (a 2D sheet), then a 2D converging section that reduces cross-sectional area smoothly can absolutely act like a throat.
Key point: it has to be a smooth contraction (think nozzle contours, not a sharp washer edge), otherwise you can get separation, turbulence, and the kind of oscillations you’ve been seeing.
In 2D, the “throat” is typically a slot throat:
throat “height” = small gap
throat “width” = large span (your flat geometry)
throat area AsubT=width×gap
At=width×gap
If it chokes, mass flow becomes primarily set by At
and upstream stagnation conditions—just like a round throat.
2) Flat expansion bells exist (conceptually): 2D C–D nozzles
The direct analogue of a conical bell is a 2D converging–diverging nozzle (sometimes called a “planar nozzle” or “rectangular nozzle” in other contexts). It expands the flow in one principal direction (height) while the width stays essentially constant.
That can work very naturally with your flat heater exit.
What you must watch:
Sidewalls: a 2D nozzle still needs side boundaries, otherwise the flow spills sideways and you lose pressure-to-axial conversion (or you generate messy 3D structures).
Boundary layer area ratio effects: a long, thin slot nozzle can be more sensitive to viscous losses than a round nozzle of the same area because the wetted perimeter per unit area can be higher.
3) “Two flat plates behind the exit” can help—but it’s not a full replacement
Two plates can provide some turning/expansion, but a real nozzle’s job is more than “something for gas to push on”:
You need the flow to expand in a controlled way (pressure gradients that accelerate primarily downstream).
If plates are too abrupt, you may get separation, shock cells, and strong unsteadiness (pogo-like behavior).
Plates can be part of a solution (as an expansion ramp / diffuser-like surface), but usually you still want a contoured ramp or contoured walls if you care about efficiency and stability.
4) Proven families of “non-circular” nozzle ideas that match your situation
Here are the ones that map best to your “flat heater → nozzle” problem:
A) 2D Slot C–D nozzle
Most direct match: keep it flat, pinch to a slot throat, then diverge with a 2D contour.
Good for: simplicity, matches your geometry, easy CFD.
Watch: side leakage, viscous losses, and making the contraction smooth.
B) Linear aerospike / expansion ramp
This is basically “half” of an aerospike wrapped into a straight line.
The expansion surface is a ramp, and the outside boundary of the jet becomes a “free boundary.”
Big advantage: less sensitive to over/under-expansion across a range of pressures.
Watch: base region behavior, shocks, and (in real life) cooling of the spike/ramp.
C) Array of small round nozzles (multi-nozzle) fed from a flat plenum
Instead of forcing the whole sheet into one circular throat, you feed many small circular throats along the width.
Advantages: uses “normal” nozzle know-how; easier manufacturing; very scalable.
Watch: manifold uniformity (each nozzle getting the same total pressure/temperature), and interactions between jets.
D) Rectangular-to-circular transition (manifold + contraction)
If your team really wants a single round nozzle, you can build a smooth 3D transition from a wide slot/plenum into a circular throat.
This is the most “classic rocket-like” compromise.
Watch: the transition must be smooth in curvature and area change, or it will behave like a loss generator (and potentially an oscillator).
5) The “is our long thin pipe acting as a throat?” idea
It can sort of behave like a restrictor, but it’s a poor substitute for a real throat because it’s dominated by friction and distributed losses rather than controlled choking and expansion.
A true throat is short and intentional: it sets
?m˙
when choked with relatively predictable behavior.A long thin pipe can:
raise upstream pressure a lot (you’re seeing that),
limit mass flow,
create conditions that encourage unsteadiness (especially if there’s compressibility + reflections + boundary-layer effects).
So: it’s not crazy that the pipe behaves “throat-like” in the sense of restricting flow—but it’s usually not the nozzle behavior you want if your goal is efficient thrust.
What I’d recommend for your next “flat-heater-first” nozzle experiment
If you want the cleanest, least controversial step:
Make a short, smooth 2D contraction right at the heater exit to a slot throat area you can compute.
Add a short 2D diverging section at your engineer’s preferred half-angle (or a gentle contour).
Put sidewalls so the flow can’t spill sideways and create weird 3D losses.
Then compare: thrust stability + mass-flow behavior + ISP.
That gives you something that is “not circular,” but still “a real nozzle” in the fluid-dynamics sense.
If you want me to tie this directly to your engineer’s chosen cone angle from the PNG, I’ll need you to re-upload that image (the earlier file handle has expired on my side).
(th)
Offline
Like button can go here
This post is about the NewNozzle CFD sprint to 16 seconds. The run completed without incident, and to my surprise, the outflow does not yet mach intake.
The system is still building up mass in the heating tube. This is pretty amazing (from my perspective).
There can't be very many researchers elsewhere on the planet studying propulsion in this specific configuration.
The only reason ** we ** are studying this is the inspired vision of kbd512, that led us to this configuration.
Here is the tail of the run to 16 seconds:
ExecutionTime = 11522.65 s ClockTime = 11606 s << 3+ hours for 8 seconds of simulation time
volFieldValue p_avg write:
volAverage(heaterCells) of p = 665639.9 << This is getting ridiculous !!! over 6 barvolFieldValue T_avg write:
volAverage(heaterCells) of T = 757.2487 << This is average over 9600 cellsvolFieldValue rho_mass write:
volIntegrate(heaterCells) of rho = 4.465979 << We are still building mass!volFieldValue T_min_global write:
min() of T = 26.83986 at location (-239.5 0.175 2.292485e-17) in cell 8880 << That 20 Kelvin gas is hitting 26+ Kelvin right at the doorvolFieldValue T_max_global write:
max() of T = 1789.563 at location (-0.5 -0.195 -2.162287e-17) in cell 239 << kbd512's goal is 3000surfaceFieldValue massFlow_in write: << We are pushing 2 kg/s in at the intake
sum("inlet") of phi = -2surfaceFieldValue massFlow_out write: << and we are NOT YET at equilibrium after 16 seconds
sum("outlet") of phi = 1.893913End
I'll stop here to make some plots ... I'm tempted to run for another 16 seconds to see if we close in on equivalence.
Something to think about is whether any material can hold together against the pressure and temperature we are seeing.
This is the specification for the geometry of NewNozzle at this moment:
Checking geometry...
Overall domain bounding box (-240 -0.2 -0.05) (5 0.2 0.05)
The measurements above are in meters, so the width is 40 centimeters, matching extendedMerlin
However, the Z dimension is slightly thinner than extendedMerlin, so 10 centimeters thick
The only thing I can think of that might account for the ** unbelievable ** thrust is the dwell time of mass in the heating tube.
If mass in the tube is 4+ kilograms, then that must mean the dwell time for each kilogram is 2 seconds, so each kilogram would be exposed to 80 MW of heating.
this is just a guess, so I'm hoping ChatGPT5 can make sense of this output.
Update: The model contains 9800 cells:
cells: 9800
Volume is 9.8 cubic meters (245*.4*.1), so each cell has a volume of "0.001" cubic meters.
I deduce that the Mesh Create utility created cells 10 centimeters on a side.
The NewNozzle Mesh has 9800 cells, and the total volume is 9.8 cubic meters.
Each cell has a volume of .001 meter or 1000 cubic centimeters.
(th)
Offline
Like button can go here
Here are the plots for the 16 second run of NewNozzle (heater without a nozzle)
The heated volume is 240 meters log, 40 centimeters wide, and 10 centimeters thick. 40 MW of heating is applied mathematically to the 9600 cells of the heater.
Hydrogen is preloaded in that volume at 1 bar and 21 Kelvin.
At start of the run, 40 MW is applied and the exit port is set to near vacuum (20,000 Pascal)
The intake ramps up to a flow of 2 kg/s over 2 seconds. Temperature at the intake is 20 Kelvin.
Thrust (unbelievable)
Pressure - does this make sense?
Mass - looks reasonable given 2 kg/s input
Temperature average over 9600 cells
Minimum temperature near intake - looks reasonable
Maximum temperature - note 5000 max cut off by safety settings
The time slices held study ... the equations were loafing.
The equations were loafing ... there was no hard work to do.
I assume there must be multiple errors impacting this run, but at this point I have no idea what they are.
(th)
Offline
Like button can go here
Once you get gas flow and it warms up a bit, you are probably "thermal choking" somewhere in the heater tube. That leads to a sonic or very low supersonic speed at the exit, which really does produce a thrust. You need either less massflow or a much larger flow cross-section area to avoid that thermal choking. It does raise the pressures inside the heater tube. Although, maybe you really want the thermal choking. I honestly dunno. -- GW
GW Johnson
McGregor, Texas
"There is nothing as expensive as a dead crew, especially one dead from a bad management decision"
Offline
Like button can go here
For GW Johnson re #270 ... thanks for taking a look at the problem!
In this post, ChatGPT5.2 is trying to figure it out .... The thrust reported is greater than physics would (apparently) support for our fluid flow and power input.
Absolutely — here is a clean BBCode version, written for a text-only forum, with no subscripts/superscripts, and structured so your engineer and team can evaluate it carefully.
Why the thrust numbers look “too good” — and what is really happening
Short answer:
We do not have a defect in OpenFOAM itself, but the current thrustCalc configuration is reporting a force that does not correspond to a real rocket thrust situation. The solver is faithfully integrating forces, but the physical interpretation is wrong for the present geometry.1. What thrustCalc is actually computing
The thrust post-processor computes force from:
[] Pressure acting on the outlet patch
[] Momentum flux leaving the outlet(Optionally) subtraction of inlet contributions, if configured
In mathematical terms (text form):
Thrust = Integral over outlet of (pressure times area normal)
plus massFlowRate times exitVelocityThis is correct only if:
[] The outlet is a real exhaust plane
[] Momentum is carried away downstreamThe vehicle is not mechanically pushing on an upstream surface
2. Your current geometry violates a key rocket assumption
Right now, your model is:
[] A long heated pipe
[] No converging throat
[] No expansion bell
[] A straight outlet into vacuumIn this configuration:
There is nothing to convert pressure into directed exhaust momentum.
The gas is being heated, expands, and pushes in both directions.
3. Where the force is actually applied
Since there is no nozzle:
[] Gas pressure builds inside the pipe
[] That pressure pushes on all solid surfacesThe only place that can react that force is the intake plate
So yes — there is a force acting on the intake plate.
But that force is:
Internal to the vehicle
It is equivalent to pushing on the inside of a closed container.
4. Why thrustCalc reports thrust anyway
thrustCalc does not know about vehicle structure.
It sees:
[] Pressure at the outlet patch
[] Net momentum flux through that patchSo it assumes:
"This must be thrust."But in your case:
[] The outlet plane is not a nozzle exit
[] The momentum is not free-stream exhaustThe pressure force is reacting internally
Result: thrustCalc reports a force that looks like thrust, but is not producing vehicle acceleration.
5. Why the magnitude looks plausible (and misleading)
Your numbers (~28 kN at 2 kg/s) appear reasonable because:
[] Heating raises pressure to ~7 bar
[] Outlet area is large (0.04 square meters)Pressure times area gives tens of kilonewtons
Numerically:
Pressure ~720,000 Pa
Area ~0.04 m2
Force ~28,800 NThe math is correct.
The physics interpretation is not.
6. Why pressure increases without a nozzle
You asked:
"Why would pressure decrease if we restrict the exit?"
Answer:
[] With no throat, flow is not choked
[] The pipe itself becomes the restriction
[] Heating adds energy faster than mass can escape
[] Pressure rises until massFlowOut approaches massFlowInThis is exactly what your logs show:
[] inlet ~2.0 kg/s
[] outlet ~1.77 kg/spressure still rising
This is expected behavior for a heated duct without a nozzle.
7. Is this producing real thrust?
No.
At present:
[] The force acts on internal surfaces
[] Momentum is not collimatedNo expansion surface exists to redirect pressure into axial thrust
This is not a rocket yet.
8. Why adding a nozzle changes everything
A nozzle does three essential things:
[] Chokes the flow at a defined throat
[] Forces pressure to convert into velocityRedirects lateral pressure into axial momentum
Only then does:
[] Pressure drop downstream
[] Exhaust carry momentum awaythrustCalc become physically meaningful
9. Conclusion for the team
[] thrustCalc is behaving correctly
[] The reported thrust is not vehicle thrust
[] The force is reacting internally on the intake plate
[] The pipe is acting as a pressure vessel, not a nozzleThe results are a diagnostic success, not a failure
Next step:
Introduce a throat (circular or slot) and an expansion surface (cone, wedge, or aerospike). Once that exists, the same thrust calculation will become physically valid.If you want, next time we can:
[] Adjust thrustCalc to explicitly subtract inlet reaction
[] Estimate the true net vehicle force before the nozzleDefine acceptance criteria for “real thrust” in your new geometry
You’ve uncovered a very important lesson here — one that many CFD users never notice.
It is possible that the explanation is in that response. I haven't had time to study it in detail.
I'm attending the monthly Linux meeting at the moment. We are looking at the alternative distributionis of Linux that are available.
(th)
Offline
Like button can go here
ChatGPT5.2 generated another version of it's reply about the no-nozzle test setup:
Thank you for investing so much of your time and resources to address this mystery. We have an opportunity to learn more about OpenFOAM and even a bit of physics as we proceed. I think your hint that the input could be involved is certainly worth investigating. I don't think you remember, but we actually discussed this very question recently. At the time, my recollection is that your instance at that time did not think it was a significant factor in whatever we were doing. This time I suspect is is a major factor. PS ... your BBCode looks terrific. Thanks for adapting to our environment. It is partly why we hacve held onto our (small) international membership all these years.
Absolutely — here is a clean, forum-ready BBCode version of the explanation, written entirely in text math, with no subscripts, no superscripts, and aligned with how thrustCalc, OpenFOAM, and physics are actually behaving in your current case.
Why the “No-Nozzle Pipe” Appears to Produce Thrust
Short answer:
No, we did not invent a new rocket nozzle — but we did uncover an important modeling artifact.
The thrust you are seeing is real force being computed correctly by OpenFOAM — but it is not jet thrust.
It is primarily pressure reaction force at the inlet.What thrustCalc is actually measuring
The thrustCalc utility computes:
[]Pressure forces on the selected patch (usually the outlet)
[]Momentum flux crossing that patchIn equation form (text only):
Thrust = integral over patch of (pressure times normal area)
plus integral over patch of (mass flow times velocity)This is correct physics — but it only tells the truth if the flow topology matches a real rocket.
What your current geometry actually is
Your current case is:
[]A long heated pipe
[]No converging throat
[]No expansion bell
[]A relatively large outlet areaMass flow imposed at the inlet
That means:
[]There is no proper jet formation
[]There is no supersonic expansionThere is no momentum collimation
Yet thrustCalc still reports ~28 kN.
Why?
Where the force really comes from
Because the exit is not choked, pressure builds in the pipe.
That pressure must be balanced somewhere.
The only solid surface that can react it is:
The inlet plate
So what OpenFOAM is modeling is essentially:
[]A pressurized pipe
[]With gas pushing on the inlet boundaryAnd that force being interpreted as vehicle thrust
This is physically equivalent to:
[]A pressure-fed piston
[]Not a rocket exhaustWhy this does NOT violate conservation of momentum
Momentum is conserved.
What happens is:
[]Mass enters with prescribed velocity
[]Pressure rises due to heating and restriction
[]That pressure exerts force on the inlet
[]The inlet reaction force is counted as thrustThis is allowed — but it is not propulsion.
A rocket requires:
[]Momentum leaving the system
[]At high directed velocityYour system currently produces force without effective exhaust momentum.
Why pressure increased instead of decreasing
You asked:
"Why would pressure decrease if the exit is restricted?"
It would not.
Correct behavior:
[]Restrict exit area
[]Pressure inside increasesMass accumulates until outlet flow matches inlet flow
This is exactly what your plots show.
The pressure rise is correct.
The thrust interpretation is not.
Why thrustCalc is not “broken”
thrustCalc is doing exactly what it was designed to do.
The problem is that:
[]The geometry does not match rocket assumptions
[]The inlet boundary is doing mechanical workThat work is being misinterpreted as thrust
This is a modeling context issue, not a bug.
Key diagnostic that confirms this
In a real rocket:
[]Most thrust comes from exit momentum
[]Pressure at inlet is irrelevantIn your case:
[]Exit velocity is modest
[]Exit momentum is too small for 28 kNPressure in heater is very high
That can only mean:
The force is coming from pressure on the inlet.
What this means for the project
This result is actually very valuable.
It proves:
[]Volume heating is working
[]Mass flow control is working
[]Pressure generation is real
[]Energy accounting is consistentBut it also proves:
[]A nozzle (or expansion surface) is absolutely required
[]Without it, thrust numbers are not meaningfulWhat changes once a nozzle is added
When you add:
[]A throat
[]An expansion wedge or bellThen:
[]Pressure converts to directed velocity
[]Momentum leaves the system
[]Inlet force becomes negligible
[]Thrust becomes physically realAt that point, thrustCalc will report something you can trust.
Bottom line
No, we did not invent a new rocket nozzle.
But we did uncover a subtle and important modeling truth:A pressurized pipe with imposed inflow can produce large forces — without being a rocket.
Your instincts were correct to question the result.
This is exactly the kind of insight that separates:
[]“The solver runs”
[]from“The physics is right”
If you want, next we can:
[]Design the expansion wedge in 2D
[]Define the equivalent 3 cm throat area correctly
[]Predict what thrust should be before running it
[]Add diagnostics to separate inlet force from exhaust momentumYou’re asking the right questions — and this is a very strong foundation for the nozzle phase.
I haven't read this yet so will have to do so later.
The Linux meeting is still in progress. We are looking at one of many distributions of Linux. This one is a descendant of Sun's Solaris.
(th)
Offline
Like button can go here
I pressed ChatGPT5.2 about the weird numbers we seeing, and it came back with a concrete proposal for how we can test the theory.
The theory is that a feature of OpenFOAM that operates automatically when certain behaviors are specified in the configuration file.
The feature in question is automatic compensation for pressure when intake is specified as 2 kg/s (or any fixed amount), In the earlier stages of learning about OpenFOAM this default behavior could be (and was) ignored. Now it is a massive contributor to the forces in play.
Yes — we absolutely can obtain the intake contribution, and you are right to focus on it.
Short answer (for the meeting):
Yes. OpenFOAM can explicitly report the force contribution from the inlet patch.
We can compute the momentum flux and pressure force acting on the inlet surface using the same mechanism already used for the outlet thrust calculation, simply applied to the inlet patch.That inlet force can be large in your current configuration (no nozzle, high internal pressure, constrained outlet), and it can fully explain the apparent “thrust” you are seeing.
When you’re ready after the meeting, the next concrete step is to:
add a second thrust/force calculation on the inlet patch, and
compare inlet force vs outlet force directly.
That comparison will tell us immediately whether thrustCalc is behaving correctly and where the force is really coming from.
Update ... the explanation that ChatGPT5.2 is considering might also account for the increasing mass in the intake pipe.
I have noticed that the output never seems to reach equilibrium with the intake. This can be be explained if the intake is pumping energy into the heating pipe automatically as pressure builds.
(th)
Offline
Like button can go here
In this post, ChatGPT5.2 addresses the question of a possible feedback loop in the current configuration of the case:
Thank you ** very ** much for the action plan!I think we can now understand why the load of mass in the heating pipe is increasing! I think we have a feedback loop in place. As the pressure inside the pipe increases, the intake automatically increases pressure in order to keep 2 kg/s coming in. What do you think?
Thought for 11sYes — that absolutely sounds like what you’ve built.
With your current inlet setup, OpenFOAM is essentially doing this:
You told it: “deliver 2 kg/s, no matter what” (via flowRateInletVelocity with a massFlowRate table).
And your inlet pressure treatment (fixedFluxPressure) + the solver will happily adjust whatever pressure gradient / inlet velocity is required to achieve that target mass flow.
So if downstream conditions make it harder to push flow out (restricted outlet, poor expansion, shocks, recirculation, “wrong nozzle,” etc.), the internal pressure rises…
…and then the inlet BC reacts by “pushing harder” to keep 2 kg/s coming in.
That is a control loop, and it can behave like positive feedback in a confined volume:
Outlet can’t pass 2 kg/s yet (or is passing less).
Mass accumulates → pressure rises in the pipe.
Inlet BC increases the driving pressure/velocity to enforce 2 kg/s anyway.
Even more mass/energy gets forced in faster than it can leave → pressure rises further.
That matches your observation: rho_mass climbing, p_avg climbing, and massFlow_out < massFlow_in for a long time.
Why I’d expect pressure to drop in a “real nozzle case”
In a healthy rocket nozzle case, the nozzle provides a clean acceleration/expansion path that:
converts chamber energy into directed momentum,
and tends to establish a stable flowfield where outflow can match inflow without needing the inlet to “fight” the domain.
In your current “heated pipe + blunt/flat outlet” setup, you don’t have that tidy expansion mechanism, so you’re closer to a pressurizing tank with a leaky outlet, not a steady nozzle.
What to do (so the CFD behaves like a real feed system)
If you want to avoid “the inlet will do anything necessary” behavior, the fix is conceptual:
Option A (most physical for an engine):
Set the inlet by total conditions (e.g., total pressure + total temperature), and let the mass flow be a result of downstream choking/conditions.
This mimics a real tank/regulator: it can only supply so much pressure.Option B (still mass-flow-driven, but prevent runaway):
Keep the mass-flow inlet, but add a cap on what the inlet is allowed to do (effectively: “don’t exceed X inlet total pressure”), or change the inlet to something that doesn’t implicitly inject unlimited “pump work” into the domain.Option C (diagnostic step):
Temporarily switch to a simpler inlet (fixed velocity or fixed pressure) and see whether the pressure runaway disappears. If it does, you’ve proven the feedback-loop hypothesis cleanly.If you paste your current 0/U inlet stanza and 0/p inlet stanza for this foamWorkroom run, I can tell you which of those options is the least disruptive edit for your exact solver/thermo combo.
I am hoping this ongoing dialog will be helpful to students who may be learning about OpenFOAM or any CFD system.
I find that there are default behaviors that help an individual to get a case running, but which contributed unwanted effects that eventually have to be recognized and addressed.
(th)
Offline
Like button can go here
A sprint by NewNozzle to 32 seconds from 16 seconds just ended after less than 13 hours.
Here are the reports at 32 seconds:
ExecutionTime = 46111.96 s ClockTime = 46555 s << less than 13 hours 46800
volFieldValue p_avg write:
volAverage(heaterCells) of p = 782376.2 << pressure still increasingvolFieldValue T_avg write:
volAverage(heaterCells) of T = 729.9388volFieldValue rho_mass write:
volIntegrate(heaterCells) of rho = 5.327563 << mass still increasingvolFieldValue T_min_global write:
min() of T = 26.83269 at location (-239.5 0.175 2.292485e-17) in cell 8880 << Temperature at intake steadyvolFieldValue T_max_global write:
max() of T = 1703.316 at location (-0.5 -0.195 -2.162287e-17) in cell 239 << Temperature at exitsurfaceFieldValue massFlow_in write: << mass flow in steady
sum("inlet") of phi = -2surfaceFieldValue massFlow_out write:
sum("outlet") of phi = 1.972983 << Still below intake after 32 seconds but closeEnd
We are approaching equivalence but "are not there yet".
Meanwhile, ChatGPT5.2 has proposed a set of data gathering and reporting functions to be installed in controlDict to look at the forces at the intake during a run. If my theory is correct, that the intake is contributing massively to the thrust reported, we should be able to subtract the intake from the total at the exit, to discover the true contribution of the heating.
It is likely we will discover that the energy required to pump hydrogen into the heater is significant, and that energy must be provided from somewhere, and the most reasonable source is from the 40 MW collected by the solar panels.
It seems likely to me that we might be able to make these changes and run again in time for next Sunday's meeting.
Here is the plot of thrust measured at the exit for the 32 second run:
(th)
Offline
Like button can go here