New Mars Forums
You are not logged in.
- Topics: Active | Unanswered
Announcement
#226 2026-01-10 14:52:29
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
This post is by itself.... It reports on ISP based upon data reported in posts 224 and 225....
In this conversation, ChatGPT5 is dealing with the noisy data that came out of the thrustCalc run that read data recorded in the postProcessing folders.
As discussed, there may be sound waves moving in that 240 meter pipe. Question: What frequency is likely? SpaceNut ... a while back, you published something here in the forum about the possible resonate column of a SpaceX Starship. kbd512's heating tube is 240 meters long. The gas inside is hydrogen, and the temperature varies inside the tube from about 20 Kelvin near the intake to ? thousands of degrees Kelvin at the throat. I have no idea how to compute the sound frequency of that pipe, but perhaps you can find a resource. You can see the peaks and values in the plot of thrust. I count 7 peaks per second. Does that match up with anything?
Here is the ISP estimate:
Yep — with the numbers you’ve provided, we can estimate ISP right now (and we can do it in a way that ignores the “flush chaos” and uses the settled portion near the end).
ISP from your 1 kg/s run (using the last ~1 second)
Definition:
Isp {s} = Thrust_N / (m_dot_kg_per_s * g0)
g0 = 9.80665 m/s^2
From your thrust_series.csv, the last ~1 second (about t = 3.05 to 4.00 s) bounces between ~9.8 kN and ~13.8 kN. Using the 20 points you showed, the mean thrust is ~11,504 N (11.5 kN).
So:
If you use m_dot = 1.000 kg/s (your intended inlet setting):
Isp ≈ 11504 / (1.000 * 9.80665) ≈ 1173 s
If you use m_dot_out = 1.020 kg/s (your tail shows massFlow_out = 1.020072):
Isp ≈ 11504 / (1.020 * 9.80665) ≈ 1150 sSo you’re in the ballpark of ~1150–1170 seconds Isp for the “late-time / quasi-steady” portion you pasted.
That’s a huge jump from ~250 s — which makes sense because with the same 40 MW and half the mass flow, you can (in principle) drive much higher exhaust velocity (until losses/limits dominate).
Why the thrust plot is “wild” in this 4-second run
Your thrust CSV has enormous spikes early (t < ~2 s), then it collapses into a much lower, noisy band (~10–14 kN). That pattern is very consistent with:
the initial inventory dump (you started with the engine full at ~1 bar and opened the outlet to vacuum), producing a strong transient,
pressure/acoustic waves bouncing through the long 240 m duct,
interaction with the porous washer / near-choked region, which can act like a “valve + resonator” and create bursts as the solution marches forward.
So: “sonic surges” (i.e., compressible-wave dynamics and/or intermittent choking/un-choking) are a very plausible explanation for the early and mid-run variability.
Practical recommendation for reporting ISP to the team
If you want a defensible number for this run without overthinking the flush:
Pick a “late-time window” where:
massFlow_in is stable at -1, and
massFlow_out is roughly steady (even if still slightly >1 due to flush).
Compute:
mean(thrust) over that window
mean(massFlow_out) over that window (or just use 1.0 if you’re explicitly evaluating the commanded inlet)
If you tell me what window you want to standardize on (e.g., last 0.5 s, last 1.0 s), I can give you a clean “reportable” ISP formula + awk one-liner consistent with your existing workflow.
(th)
Offline
Like button can go here
#227 2026-01-10 15:05:33
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
I asked ChatGPT5 to ** think ** about acoustics in our running model. I don't understand 90% of this reply, and that's being optimistic.
However someone in the membership may understand this discussion and be able to provide advice on what to do to see the sound waves that (presumably) are present. I gather that the ParaView program cannot show sound waves by default, but there may be some way to collect data that can be translated into visibility of compression waves in the medium inside the heating tube.
1. What “sonic waves” really are in your simulation
In OpenFOAM terms, what you are seeing is unsteady compressible behavior, which can include:
Compression waves (pressure waves)
Expansion waves
Shocklets / weak shocks
Choked-flow oscillations
Acoustic ringing between the washer and downstream expansion
These will show up numerically as:
Rapid oscillations in thrust
Fluctuations in pressure
Variations in massFlow_out
Local spikes in velocity magnitude
High gradients in Mach number
Your thrust CSV strongly suggests this:
Rapid jumps (10 kN → 70 kN → 10 kN)
Non-smooth oscillatory behavior
Strong sensitivity during early seconds
That is classic compressible transient behavior, not numerical noise.
2. Will ParaView show them automatically?
No.
By default, ParaView hides sonic phenomena unless you explicitly reveal them.If you only look at:
Surface
Solid Color
Raw U magnitude
…you will miss most wave behavior.
3. The best ParaView views for sonic waves (ranked)
✅ #1: Mach number (MOST IMPORTANT)If you only do one thing, do this.
How:
Apply Calculator
Create a new field:
Mach = mag(U) / sqrt(gamma * R * T)
For hydrogen (approx):
gamma ≈ 1.4
R ≈ 4124 J/(kg·K)
Color by Mach
Use Surface or Slice
What you will see:
Choked regions (Mach ≈ 1)
Supersonic pockets (Mach > 1)
Moving Mach fronts = sonic waves
This is the diagnostic for your washer test.
✅ #2: Pressure gradients (shock visualization)
How:Apply Gradient Of Unstructured Data Set
Input: p
Visualize |∇p|
What you’ll see:
Thin bright bands = compression fronts
Moving structures = pressure waves
Shocks don’t always show up as sharp “walls” in hydrogen at these conditions—often they appear as traveling pressure ridges.
✅ #3: Surface LIC (you already found gold here)
You already discovered this one instinctively—and that’s impressive.
Surface LIC shows direction and coherence of flow.
Sonic phenomena appear as:
Flow-line bunching
Rapid directional flipping
Standing-wave patterns downstream of the washer
If LIC “shimmers” or reorganizes rapidly in time → acoustic activity is present.
✅ #4: Time animation (essential)
Static images are misleading.
Do this:
Animate pressure, Mach, or |U|
Play at real time or slowed down
Watch for:
Pulses originating at the washer
Reflections off the nozzle
Standing oscillations
If sonic waves are present, you will see motion, not just gradients.
4. Why your thrust plot looks “wild” at 1 kg/s
This is the key physical explanation:
At 1 kg/s with 40 MW:
Energy per unit mass doubles
Gas heats faster
Local sound speed rises sharply
Flow becomes highly compressible
Choking + unsteady expansion is likely
Your washer creates:
A control volume
A pressure reservoir
A reflective boundary for waves
That combination is not steady by default.
What you are seeing is very likely:
Transient choking
Acoustic ringing
Pressure wave cycling
This is real physics, not a bug.
5. Why you did NOT see this at 2 kg/s
At 2 kg/s:
Lower specific energy
Lower temperature
Higher damping
Flow is “heavier” and smoother
At 1 kg/s:
You crossed into a different regime
This is exactly what you wanted to explore.
6. Should you worry?
No.
You should be encouraged.You have reached a regime where:
ISP is increasing
Compressible physics dominates
Flow control becomes design-critical
This is where:
Nozzle shaping
Washer geometry
Plenum length
Throat treatment
…all start to matter.
That means your model is doing its job.
7. One simple test to confirm sonic waves (high confidence)
Plot Mach vs time at a single probe point just downstream of the washer.
If Mach oscillates around ~1 → you have unsteady choking.
Final takeaway
Yes—ParaView can show sonic waves, and your model almost certainly contains them.
You just need to look at Mach, pressure gradients, and time animation, not static velocity plots.You are no longer debugging.
You are doing propulsion physics.When you’re ready, next logical steps would be:
Quantify Mach oscillation frequency
Correlate with thrust oscillations
Decide whether to damp or exploit the behavior
That’s an exciting place to be.
(th)
Offline
Like button can go here
#228 2026-01-10 15:51:59
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
This post is a follow up on my conversation with SpaceNut ... I asked ChatGPT5.2 to **think** about the resonant frequency we are seeing of about 7 Hertz. It found that this frequency might be a second harmonic or a first harmonic, depending upon heating conditions in the pipe.
What I observed while running ParaView (animation viewer) is that the pulses of energy are very dramatic when viewed in animation at the engine compartment. The flume from the nozzle expands and contracts in rhythm with the bursts of energy visible in the plot.
Thank you for your generous response to our question about sound waves active in our model. I've posted your explanation and hope that someone follows up with a suggestion or a request. It seems to me this is not a desirable situation and that we might be able to do something. Can you estimate the frequency of the sound? We can see 7 peaks a second. Does that correspond to the 240 meter pipe if hydrogen is the gas?
Thought for 15sYes — 7 peaks/second (~7 Hz) is absolutely in the right ballpark for an acoustic (standing-wave) mode in a 240 m hydrogen duct, depending on:
speed of sound in the gas (which depends strongly on temperature), and
whether the pipe behaves more like open–open or open–closed (and what the “effective acoustic length” really is).
1) Speed of sound in hot hydrogen (order-of-magnitude)
For a perfect gas, the sound speed is:
a = sqrt(gamma * R * T)
For hydrogen:
R (specific) ≈ 4124 J/(kg·K)
gamma ≈ 1.4 (reasonable first pass for diatomic gas; it can drift with temperature/dissociation)
So:
At T = 700 K:
a ≈ sqrt(1.4 * 4124 * 700) ≈ ~2000 m/sAt T = 1400 K:
a ≈ ~2850 m/sYour logs show T_avg ~ 700–1400 K, and T_max can be much higher, so 2000–3000 m/s is a very plausible range.
2) What frequency does a 240 m “organ pipe” give?
Let L = 240 m.
Open–open fundamental
f1 = a / (2L)
If a = 2000 m/s → f1 ≈ 2000 / 480 ≈ 4.17 Hz
If a = 2850 m/s → f1 ≈ 2850 / 480 ≈ 5.94 Hz
Open–closed fundamental
f1 = a / (4L)
If a = 2000 m/s → 2.08 Hz
If a = 2850 m/s → 2.97 Hz
So 7 Hz is:
a bit higher than open–open fundamental for ~1400 K,
but very consistent with either:
a hotter effective temperature, or
a shorter effective acoustic length, or
a higher mode (harmonic) rather than the fundamental.
3) If your measured frequency is 7 Hz, what does it imply?
If it’s open–open-like, then:
a ≈ 2 L f = 2 * 240 * 7 = 3360 m/s
That corresponds to an “effective temperature” of roughly:
T ≈ a^2 / (gamma * R)
T ≈ 3360^2 / (1.4 * 4124)
T ≈ 11,289,600 / 5773.6 ≈ ~1950 KThat’s not crazy at all given you’re seeing hot regions (high T_max) and strong gradients. In other words: 7 Hz can match a 240 m duct if the oscillation “lives” in the hotter part of the system (or the mode is not spanning the entire 240 m uniformly).
4) Why 7 peaks/sec might not be “the 240 m fundamental”
A few very plausible reasons:
Effective length is not 240 m. The acoustic cavity may be from the washer-to-nozzle region, or some subset bounded by a reflecting feature (washer restriction, area change, nozzle throat, etc.).
It’s a higher harmonic. If the fundamental is ~3.5 Hz, the 2nd harmonic is ~7 Hz. That happens easily when boundary reflections are strong.
It’s a coupled “engine” mode, not purely the straight pipe (nozzle + chamber + restriction can form a more complex resonator).
Numerics can imprint periodicity. If shocks form and move, or if waveTransmissive/outlet behavior partially reflects, you can get a “ringing” that looks acoustic even when it’s really a flow-acoustic + boundary-condition interaction.
5) How to confirm it’s acoustic (quick, decisive checks)
If you want to prove it’s an acoustic mode:
Pick 1–3 probe points (pressure is best), record p(t) at high cadence.
Do an FFT of p(t).
Compare the dominant peak to a/(2L) or a/(4L) using a reasonable T for the region where the wave lives.
If you see a strong spectral line at ~7 Hz (and maybe harmonics), it’s very likely a resonance.
Bottom line
Yes: seeing ~7 Hz is consistent with a resonance involving a ~240 m hydrogen system if the effective sound speed is in the ~3000 m/s range (hot gas) or the mode is a harmonic / shorter effective cavity. Your “7 peaks per second” observation is the right kind of number for this scale.
If you tell me whether your outlet behaves closer to “open” acoustically (low reflection) or is reflecting strongly (often true in practice), I can tell you which boundary model (open–open vs open–closed) is the better first approximation and what harmonic structure to expect.
(th)
Offline
Like button can go here
#229 2026-01-11 10:40:36
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
This post includes links to images showing the extendedMerlin engine compartment at 4 seconds of a run using 1 kg/s of hydrogen and 40 MW of energy.
The fixed images do not show the (to me interesting) "pogo" pulsing at about 7 Hertz that can be seen in the animation.
The pulses are clearly visible in the thrust plot included in a post above.
Pressure
Temperature
Velocity (pulses are most pronounced in this mode)
(th)
Offline
Like button can go here
#230 2026-01-11 22:11:35
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
Following this evening Google Meeting, during which GW Johnson evaluated results of the 1 kg/s run, I decided to begin next week's sessions with a run using 2 kg/s and no other changes. This run will last for at least 3 days, so it will be a while before we see the thrust plot. however, I'll try to catch the system at 1 second intervals to generate plots showing how things are going, as we did with the 1 kg/s case. Here is the opening status report:
ExecutionTime = 86.19885 s ClockTime = 87 s
volFieldValue p_avg write:
volAverage(heaterCells) of p = 106571.6 << We are starting at 1 bar with 14 kg hydrogen filling the modelvolFieldValue T_avg write:
volAverage(heaterCells) of T = 22.2374 << The starting temperature is 21 Kelving but heat is already pouring involFieldValue rho_mass write:
volIntegrate(heaterCells) of rho = 14.06393 << The model starts with 14 kg of hydrogen at 1 bar and 21 KelvinvolFieldValue T_min_global write:
min() of T = 19.86342 at location (4.157199 -0.4660463 0) in cell 10521 << That low temperature is near exit to vacuumvolFieldValue T_max_global write: << High temperature is near the intake
max() of T = 22.23846 at location (-240.025 -0.192 0) in cell 37549<< The location of high temperature will move to the right and end up just ahead of the throat
Courant Number mean: 0.000744156 max: 0.01669082 <<
deltaT = 1e-05
Time = 0.0063575s
In this run it will be interesting to see if oscillation occurs as we observed with 1 kg/s. In this run we are using the same flat washer for a throat, and that means we are NOT providing smooth flow for the gas, so turbulence is likely. We will have these two runs available for comparison with runs using the new mesh.
That new mesh will be difficult to integrate into the existing model. It won't necessarily be difficult to create the model of the nozzle GW Johnson is recommending, because the components appear to be shapes that are readily available in Blender (or any 3D package). What I do not know how to do is to impose the new mesh upon the existing mesh.
(th)
Offline
Like button can go here
#231 2026-01-12 16:06:42
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
This post contains links to images showing plots of performance of the extendedMerlin model at 2 seconds into a 4 second run.
Conditions for this run are identical to the previous run with 1 kg/s of feed, with the sole exception we returned to 2 kg/s ramping up over 2 seconds.
In these plots it is possible to see the initial load of 14 kg of hydrogen flushing while the intake ramps up.
A curious artifact of the plots is that temperature shows a vertical rise in the last frame. This may be caused by the procedure to generate the plots. Data for the plots was collected while the run was in progress, so it is possible that incorrect data was collected in the last frame.
Here are the 8 plots, omitting one that only fills at the end of the run:
Time between computations falls as complexity increases
Equations work harder and settle at the run setting of .7
Temperature of the heat pipe increases with 40 MW of heating
Mass decreases from 14 initial load to 2 kg at the intake
Pressure rises with flush of initial load and then falls toward 2 kg/s level
Minimum temperature falls and then settles at 20 Kelvin intake temperature
Maximum temperature climbs as flush ends. The last frame may be an artifact
(th)
Offline
Like button can go here
#232 2026-01-14 08:51:23
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
We have reached the planned stopping point at four seconds, using 2 kg/s feed with 40 MW of power and the washer.
We see the same pulsing that we saw with the 1 kg/s feed. Since the only change from the previous tests is use of the washer, I suspect the inappropriate shape of the washer is causing the pulsing.
Here is the plot for four seconds:
Here is the End of Run report:
ExecutionTime = 196472.5 s ClockTime = 196780 s << 56 hours +
volFieldValue p_avg write:
volAverage(heaterCells) of p = 329640.2 << just over 3 barvolFieldValue T_avg write:
volAverage(heaterCells) of T = 749.6522volFieldValue rho_mass write:
volIntegrate(heaterCells) of rho = 2.911246 << Nearly 3 kg in the heatervolFieldValue T_min_global write:
min() of T = 22.45924 at location (-240.275 -0.16 0) in cell 37621 << As expected with intake at 20 KelvinvolFieldValue T_max_global write:
max() of T = 1637.984 at location (-0.182823 -0.1878 3.407286e-18) in cell 3401 << Maximum just before throatsurfaceFieldValue massFlow_in write: << intake is 2 kg/s as planned
sum("inlet") of phi = -2surfaceFieldValue massFlow_out write: << Outflow should be 2 kg/s on average.
sum("outlet") of phi = 1.748118 << However, because of pulsing, it must be varyingEnd
My guess is that the outlet flow measurement was taken at an instant so we get the reading at that moment.
(th)
Offline
Like button can go here
#233 2026-01-14 09:30:19
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
This post will show the pressure plot of the 1 kg/s run vs the 2 kg/s run.
1 kg/s run to 4 seconds:
2 kg/s run to 4 seconds
The greater pressure is due to increased mass flow.
Nothing else changed.
(th)
Offline
Like button can go here
#234 2026-01-14 09:54:39
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
This post will show the mass accumulation for the 1 kg/s run vs the 2 kg/s run
1 kg/s run
Mass in system stays below 2 kg. At 1 kg/s of intake, accumulation apparently does not occur.
2 kg/s run:
Mass in system rises above 2 kg. With 2 kg/s in intake, mass accumulates due to restriction of the washer acting at the throat.
I am pretty sure this plot shows mass sticking around in the heater. Outflow should match intake but that plot did not work.
Offline
Like button can go here
#235 2026-01-14 10:15:54
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
This post contains links to other plots from the 2 kg/s run to 4 seconds:
The rise in maximum temperature may coincide with the pulses in thrust
Average temperature in the 24,000 cells of the heat pipe
Minimum temperature is set by intake of 20 Kelvin
Possible coincidence with pulsing
Equations are limited by run parameter of 7 so time decreases to allow work to be completed
All cells in the model are processed at each time interval
Multiple equations are solved for each cell, and many equations are executed multiple times.
(th)
Offline
Like button can go here
#236 2026-01-15 19:10:13
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
This post is the first of what I hope will be a series about creating an OpenFOAM mesh from scratch, using OpenFOAM tools.
For the past 8 months we have been working with an OpenFOAM model of the SpaceX Merlin engine. I don't know anything about the origin of the model, but in the past few days I've come around to the view that this may have been a student project. I trust the student got a good grade for the work, because it has withstood major modification for the Optical Plane application, and after months of learning how to run the case, it is purring like a kitten. Even the indignity of forcing a flat washer into the throat did not break it, although that action ** did ** cause it to show tremors instead of the nice steady flow with the original smooth throat. We will continue to study the inherited model, in hopes of learning how to remove some of the features that appear to be left over from the student tutorial that may have been used to build the Merlin engine.
In ** this ** post I will show a file that (I am told) is designed to begin the process of building the new mesh. For those students who might run across this topic while scanning the Internet, we have created a directory structure that looks like this:
Top level of case where runs are requested: foamWorkroom/
Subdirectories at this level: system, constant and 0 << Zero ... as in Zero time. During a run OpenFOAM creates time directories full of data
The constant directory only has one subdirectory at this point: polyMesh << Empty at this point.
My advisor in this learning exercise is ChatGPT5.2. We (and ChatGPT4o) have been working without interruption for 8+ months. Due to limitations of how much memory OpenAI can allocate to individual customers, ChatGPT5.2 cannot remember everything we have done, and as a human I forget many details. However, we've done a decent job of saving transcripts. My main contribution is to manage the flow of work. During the early phase,we lost a lot of time because I did not understand how to manage a co-worker of this type. Things are MUCH more stable these days. I handle ** all ** updates to the system, and generally I only allow one change at a time. For example, we are currently running the extendedMerlin model with 2 kg/s of intake at 40 MW of power, from 4 seconds to 8 seconds. The purpose of this sprint is to evaluate the performance of the model after most of the startup fluctuations are over with. So, while ** that ** going on, I'm starting to work on the new mesh.
The name of this file is: blockMeshDict, and it is located in the system folder.
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / Field | OpenFOAM: The Open Source CFD Toolbox |
| \\ / Operation | Version: v12 |
| \\ / And | |
| \\/ Manipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
convertToMeters 1;
// ---- Geometry knobs (edit these) ----
xmin -245; // start of inlet
xsplit 0; // split between "heatedWall" and "wallSection"
xmax 5; // outlet plane
ymin -0.20;
ymax 0.20;
zmin -0.05;
zmax 0.05;
// ---- Mesh resolution knobs (edit these) ----
// Keep nz=1 for 2D empty. Increase nx1/nx2 and ny as needed.
nx1 240; // cells in heater region (xmin->xsplit)
nx2 5; // cells in short downstream region (xsplit->xmax)
ny 40;
nz 1;
vertices
(
// Block 1: xmin -> xsplit
( $xmin $ymin $zmin ) // 0
( $xsplit $ymin $zmin ) // 1
( $xsplit $ymax $zmin ) // 2
( $xmin $ymax $zmin ) // 3
( $xmin $ymin $zmax ) // 4
( $xsplit $ymin $zmax ) // 5
( $xsplit $ymax $zmax ) // 6
( $xmin $ymax $zmax ) // 7
// Block 2: xsplit -> xmax
( $xmax $ymin $zmin ) // 8
( $xmax $ymax $zmin ) // 9
( $xmax $ymin $zmax ) // 10
( $xmax $ymax $zmax ) // 11
);
blocks
(
hex (0 1 2 3 4 5 6 7) ( $nx1 $ny $nz ) simpleGrading (1 1 1)
hex (1 8 9 2 5 10 11 6) ( $nx2 $ny $nz ) simpleGrading (1 1 1)
);
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(0 3 7 4)
);
}
outlet
{
type patch;
faces
(
(8 10 11 9)
);
}
// “Heated” walls on Block 1 (top + bottom)
heatedWall
{
type wall;
faces
(
// bottom wall of block 1
(0 4 5 1)
// top wall of block 1
(3 2 6 7)
);
}
// “Non-heated” walls on Block 2 (top + bottom)
wallSection
{
type wall;
faces
(
// bottom wall of block 2
(1 5 10 8)
// top wall of block 2
(2 9 11 6)
);
}
// 2D empty planes (z-min and z-max) across both blocks
frontAndBackPlanes
{
type empty;
faces
(
// z-min faces
(0 1 2 3)
(1 8 9 2)
// z-max faces
(4 7 6 5)
(5 6 11 10)
);
}
);
mergePatchPairs
(
);Please note that the new mesh is going to use the same component names (eg, outlet) as we are using in the current running model.
The reason is (as I understand it) that OpenFOAM only cares about the named components. They must be present.
Components can have any name (I gather) but since the mesh we are replacing has these component names, they must be present in the new model.
Please note as well that we have chosen to deploy the new mesh with the intake heater pipe as the anchor for the new model. We've extended the heater by 5 meters to account for the length of the engine compartment in the current running model. Later, when the model is working with just the heating pipe, we'll add the new nozzle. The purpose of this sequence is to make sure the existing case is comfortable with the new minimal mesh. If all goes to plan, 2 kg/s of hydrogen will enter the intake as it does now. The heating subsystem will heat the gas as it moves through the pipe. At the end of the pipe, the gas will simply gush out into vacuum. I would not expect much pressure to build up but it might due to heating.
In the next post, I plan to show the transcript as I attempt to run this new command file.
(th)
Offline
Like button can go here
#237 2026-01-15 19:27:21
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
Assuming a student finds this series on OpenFOAM, the failures should be as valuable as the successes. In this case, we stumbled right out of the gate, because a control file called "controlDict" is not present. That file is not present because we are starting the case from scratch, with empty directories.
blockMesh
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
Build : 12-86e126a7bc4d
Exec : blockMesh
Date : Jan 15 2026
Time : 20:12:43
Host : "tom-GA-78LMT-01"
PID : 2293233
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 time--> FOAM FATAL ERROR:
cannot find file "/home/tahanson/OpenFOAM/tahanson-v2306/run/foamWorkroom/system/controlDict"From function virtual Foam::autoPtr<Foam::ISstream> Foam::fileOperations::uncollatedFileOperation::readStream(Foam::regIOobject&, const Foam::fileName&, const Foam::word&, bool) const
in file global/fileOperations/uncollatedFileOperation/uncollatedFileOperation.C at line 539.FOAM exiting
I've notified ChatGPT5.2 of the situation, and am waiting for a response.
My experience with ChatGPT over the months has been that responsiveness varies greatly. My experience right now is that it can take half an hour or more for ChatGPT5.2 to turn around my request, and thinking time is frequently measured in tens of seconds. I assume this must be in part due to the length of time we've been working on this project, but I expect it must also be due in part to the heavy workload of customers around the planet.
This is a work day, and college students are back in school, so heavy usage is to be expected.
I'm not in a super hurry at this point. we have a working model and it is running.
This particular run will take 72 hours minimum, and the latest estimate indicated we can expect at least 4 hours beyond that initial estimate.
The reason for the time is the large number of cells in the model (about 37000) and the difficulty of computing the state of each cell every time slice. Time slices are shrunk as the work gets more challenging. Work gets more challenging with heating.. Running the same flow at 20 Kelvin resulted in 4 hours per second of simulation time. We are now seeing 20~ hours of wall clock time for every second of sim time.
(th)
Offline
Like button can go here
#238 2026-01-15 19:57:32
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
Since ChatGPT5.2 is busy, I asked Gemini for advice, and it suggested several choices, but the one I like was to simply copy our existing controlDict. I did that and blockMesh ran. Here is the output.
blockMesh
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
Build : 12-86e126a7bc4d
Exec : blockMesh
Date : Jan 15 2026
Time : 20:51:40
Host : "tom-GA-78LMT-01"
PID : 2314152
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 timeDeleting polyMesh directory
"/home/tahanson/OpenFOAM/tahanson-v2306/run/foamWorkroom/constant/polyMesh"
Reading "blockMeshDict"Creating block mesh from
"system/blockMeshDict"
Creating block edges
No non-planar block faces defined
Creating topology blocks
Creating topology patchesCreating block mesh topology
Check topology
Basic statistics
Number of internal faces : 1
Number of boundary faces : 10
Number of defined boundary faces : 10
Number of undefined boundary faces : 0
Checking patch -> block consistencyCreating block offsets
Creating merge list .Creating polyMesh from blockMesh
Creating patches
Creating cells
Creating points with scale 1
Block 0 cell size :
i : 1.020833 .. 1.020833
j : 0.01 .. 0.01
k : 0.1
Block 1 cell size :
i : 1
j : 0.01 .. 0.01
k : 0.1Writing polyMesh
----------------
Mesh Information
----------------
boundingBox: (-245 -0.2 -0.05) (5 0.2 0.05)
nPoints: 20172
nCells: 9800
nFaces: 39485
nInternalFaces: 19315
----------------
Patches
----------------
patch 0 (start: 19315 size: 40) name: inlet
patch 1 (start: 19355 size: 40) name: outlet
patch 2 (start: 19395 size: 480) name: heatedWall
patch 3 (start: 19875 size: 10) name: wallSection
patch 4 (start: 19885 size: 19600) name: frontAndBackPlanesEnd
(th)
Offline
Like button can go here
#239 2026-01-15 20:23:41
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
checkMesh is an OpenFOAM utility that will examine our new mesh and report on it's condition.
checkMesh
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
Build : 12-86e126a7bc4d
Exec : checkMesh
Date : Jan 15 2026
Time : 21:07:31
Host : "tom-GA-78LMT-01"
PID : 2314397
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 polyMesh for time = 0
Time = 0s
Mesh stats
points: 20172
internal points: 0
faces: 39485
internal faces: 19315
cells: 9800
faces per cell: 6
boundary patches: 5
point zones: 0
face zones: 0
cell zones: 0Overall number of cells of each type:
hexahedra: 9800
prisms: 0
wedges: 0
pyramids: 0
tet wedges: 0
tetrahedra: 0
polyhedra: 0Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology
inlet 40 82 ok (non-closed singly connected)
outlet 40 82 ok (non-closed singly connected)
heatedWall 480 964 ok (non-closed singly connected)
wallSection 10 24 ok (non-closed singly connected)
frontAndBackPlanes 19600 20172 ok (non-closed singly connected)Checking geometry...
Overall domain bounding box (-245 -0.2 -0.05) (5 0.2 0.05)
Mesh has 2 geometric (non-empty/wedge) directions (1 1 0)
Mesh has 2 solution (non-empty) directions (1 1 0)
All edges aligned with or perpendicular to non-empty directions.
Max cell openness = 3.252607e-16 OK.
Max aspect ratio = 102.0833 OK.
Minimum face area = 0.001. Maximum face area = 0.1020833. Face area magnitudes OK.
Min volume = 0.001. Max volume = 0.001020833. Total volume = 10. Cell volumes OK.
Mesh non-orthogonality Max: 0 average: 0
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 2.227334e-13 OK.
Coupled point location match (average 0) OK.Mesh OK.
End
It appears the startup mesh is in good condition. However, it is too long by 5 meters so I'll edit the blockMesh and re-run it.
(th)
Offline
Like button can go here
#240 2026-01-15 20:42:55
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
Since ChatGPT5.2 is bogged down with customers, I'm happy to report that Gemini stepped in to pick up the slack. It seems to understand what needs to be done to add functionality to the new case. The mesh itself appears to be in excellent condition. i did rerun blockMesh to correct the length of the intake heater from 245 to 240, and that procedure appears to have gone well.
In the next phase of the transition, I ** think ** that we are going to have to visit each of the control files to adjust them for the new mesh.
Update Friday:
The current run of extendedMerlin is proceeding well. It just crossed 7 seconds and will finish around noon on Saturday, at 8 seconds. I am looking forward to seeing the results. We are feeding 2 kg/s and providing 40 MW of heating while impeding exit modestly with a flat washer shaped fixture in the throat. We expect to see pulsations due to the incorrect shape of the throat. What I am more interested in is the ISP we achieve with this configuration.
The new mesh development is under way. In the next post I will begin a series that I hope will be of value to someone studying OpenFOAM. My plan is to report configuration of a case for the new mesh, in a step-by-step manner with the actual files and discussion of the contents of the files.
Here is the status of the run at 7 seconds:
ExecutionTime = 37269.09 s ClockTime = 37291 s
volFieldValue p_avg write:
volAverage(heaterCells) of p = 334755.6 << Over 3 bar in the heating pipevolFieldValue T_avg write:
volAverage(heaterCells) of T = 716.9687 << Average temperature in 24,000 cellsvolFieldValue rho_mass write:
volIntegrate(heaterCells) of rho = 3.008885 << Mass in the heating pipe. 2 kg feed + 1 kg backed up by throatvolFieldValue T_min_global write:
min() of T = 22.45848 at location (-240.275 -0.16 0) in cell 37621 << Intake temperature is 20 KelvinvolFieldValue T_max_global write:
max() of T = 1545.734 at location (-0.182823 -0.1878 3.407286e-18) in cell 3401 << Hottest at just before throatCourant Number mean: 0.02224844 max: 0.6999971
deltaT = 2.085652e-06
Time = 7.0389087s << 7 seconds into 8 second run
(th)
Offline
Like button can go here
#241 2026-01-16 20:44:23
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
The purpose of the next series of posts is to provide a tutorial for someone who might be studying OpenFOAM.
We are working with a fully debugged configuration that was created to study the original Merlin engine, and which was adapted to develop a design that might be considered for an Optical Plane deep space vessel. We have a topic about that initiative, if anyone is interested.
Our current initiative is to replace the mesh with one that is better suited for hydrogen. The current throat diameter is 25 centimeters, and GW Johnson has recommended 3 centimeters as throat diameter for the next model. kbd512 found this diameter in numerous existing rockets in the open market, with numerous working examples deployed in space.
In prior posts we have shown how the new mesh was created, and what the checkMesh looks like. Now we are building the case.
A case has a fixed structure. Here I will paste part of a tree of the existing, running case.
├── 7.0500007 << Time folder created by current run
│ ├── Co << The solver computes data about each cell in the model at each time interval
│ ├── p << The write interval in this run is every .05 second of simulation time
│ ├── phi << The number of computation intervals between write events varies
│ ├── rho << The length of a computation time interval decreases with complexity of the situation
│ ├── T << In general, the hotter the gas in the model, the longer a run takes
│ ├── U
│ └── uniform
│ └── time
├── applications << This folder contains software added to the model to compute thrust at the end of a run
│ └── utilities
│ └── thrustCalc
│ ├── Make
│ │ ├── files
│ │ ├── linux64GccDPInt32Opt
│ │ │ ├── files
│ │ │ ├── thrustCalc.C.dep
│ │ │ └── thrustCalc.o
│ │ └── options
│ └── thrustCalc.C
├── constant << This is a required folder for any OpenFOAM case
│ ├── fvModels << This control file allows use of "volume heating"
│ ├── physicalProperties << This control file covers gas or fluid used in the case
│ ├── polyMesh << This folder inside constant holds the mesh and it's associated control files
│ │ ├── boundary
│ │ ├── cellZones
│ │ ├── faces
│ │ ├── faceZones
│ │ ├── neighbour
│ │ ├── owner
│ │ ├── points
│ │ └── sets
│ │ ├── addedCells
│ │ ├── ALL_heated
│ │ ├── heatedWallFaces
│ │ ├── heaterCells
│ │ ├── heat_fromCells
│ │ ├── HEAT_minusY
│ │ ├── HEAT_plusY
│ │ ├── HEAT_set
│ │ ├── inletFaces
│ │ ├── inletPick
│ │ ├── inletZone
│ │ ├── oBotPick
│ │ ├── oTopPick
│ │ ├── outletPick
│ │ ├── wallFaces
│ │ ├── washerCells
│ │ └── washerZone
│ ├── thermophysicalProperties << This is a required control file
│ └── turbulenceProperties << This is an optional control file if turbulence is expected
├── howtonotes
│ ├── CenterZoomScript.py
│ ├── historyAsOf20250922
│ ├── HowToCheckMass-for-f-in-postProcessing.txt
│ ├── HowToNotesRightShift.txt
│ ├── HowToRemove0directories
│ ├── HowToTellIfRunStopped.txt
│ ├── savePrint-controlDict.txt
│ └── ShowEngineCompartment.py
├── images
├── log72.Washer-2kg-0-4 << This run started with a sprint from 0 to 4 seconds
├── log73.Washer-2kg-4-8 << We started another run but interrupted the run to 8 seconds
├── log74.Washer-2kg-5-8 << This restart ended early due to a system error.
├── log75.Washer-2kg-6-8 << This restart is in progress, with planned end point of 8 seconds
├── postProcessing << This folder is created by OpenFOAM. It contains data that can be examined after a run
│ ├── heatBoxStats
│ │ ├── 0
│ │ │ └── volFieldValue.dat
│ │ ├── 4
│ │ │ └── volFieldValue.dat
│ │ ├── 5.75
│ │ │ └── volFieldValue.dat
│ │ └── 6.45
│ │ └── volFieldValue.dat
│ ├── massFlow_in
│ │ ├── 0
│ │ │ └── surfaceFieldValue.dat
│ │ ├── 4
│ │ │ └── surfaceFieldValue.dat
│ │ ├── 5.75
│ │ │ └── surfaceFieldValue.dat
│ │ └── 6.45
│ │ └── surfaceFieldValue.dat
│ ├── massFlow_out
│ │ ├── 0
│ │ │ └── surfaceFieldValue.dat
│ │ ├── 4
│ │ │ └── surfaceFieldValue.dat
│ │ ├── 5.75
│ │ │ └── surfaceFieldValue.dat
│ │ └── 6.45
│ │ └── surfaceFieldValue.dat
│ ├── p_avg
│ │ ├── 0
│ │ │ └── volFieldValue.dat
│ │ ├── 4
│ │ │ └── volFieldValue.dat
│ │ ├── 5.75
│ │ │ └── volFieldValue.dat
│ │ └── 6.45
│ │ └── volFieldValue.dat
│ ├── rho_mass
│ │ ├── 0
│ │ │ └── volFieldValue.dat
│ │ ├── 4
│ │ │ └── volFieldValue.dat
│ │ ├── 5.75
│ │ │ └── volFieldValue.dat
│ │ └── 6.45
│ │ └── volFieldValue.dat
│ ├── T_avg
│ │ ├── 0
│ │ │ └── volFieldValue.dat
│ │ ├── 4
│ │ │ └── volFieldValue.dat
│ │ ├── 5.75
│ │ │ └── volFieldValue.dat
│ │ └── 6.45
│ │ └── volFieldValue.dat
│ ├── thrustCalc
│ │ └── outlet
│ ├── T_max_global
│ │ ├── 0
│ │ │ └── volFieldValue.dat
│ │ ├── 4
│ │ │ └── volFieldValue.dat
│ │ ├── 5.75
│ │ │ └── volFieldValue.dat
│ │ └── 6.45
│ │ └── volFieldValue.dat
│ └── T_min_global
│ ├── 0
│ │ └── volFieldValue.dat
│ ├── 4
│ │ └── volFieldValue.dat
│ ├── 5.75
│ │ └── volFieldValue.dat
│ └── 6.45
│ └── volFieldValue.dat
├── system << This is a required folder in an OpenFOAM case
│ ├── controlDict << This is a required file ... it is used to start and control a run
│ ├── createPatchDict
│ ├── extrudeMeshDict
│ ├── fvConstraints << This is a required file. This file contains parameters that cover edge cases to a run can proceed
│ ├── fvSchemes
│ ├── fvSolution
│ ├── renumberMeshDict
│ ├── setFieldsDict
│ ├── thrustCalcDict
│ ├── topoSetDict
│ └── topoSetDict.washer
└── << The list above contains many files that were created over 8 months of fine tuning339 directories, 1084 files << This is a running case, so time folders were not included in the report
<< Time folder 0 (zero) is required for every OpenFOAM case
Here is the tree of the new case as it stands at this moment:
├── 0 << This is a required folder in an OpenFOAM case
│ ├── Co
│ ├── p << This is a required file in an OpenFOAM case
│ ├── T << This is a required file in an OpenFOAM case
│ ├── U << This is a required file in an OpenFOAM case
│ ├── uniform
│ │ └── time
│ └── wallHeatFlux
├── constant << This is a required folder in an OpenFOAM case
│ └── polyMesh << This is a required folder in an OpenFOAM case
│ ├── boundary << This is a generated file in an OpenFOAM case
│ ├── faces << This is a generated file in an OpenFOAM case
│ ├── neighbour << This is a generated file in an OpenFOAM case
│ ├── owner << This is a generated file in an OpenFOAM case
│ └── points << This is a generated file in an OpenFOAM case
└── system << This is a required folder in an OpenFOAM case
├── blockMeshDict << This control file was configured to create the new mesh
└── controlDict << This is a required file in an OpenFOAM case6 directories, 13 files
We are adding files one by one or in small groups. We are currently working with the 0 (time) folder
Posts following this will will contain the actual files and text discussing them
(th)
Offline
Like button can go here
#242 2026-01-16 21:31:53
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
After we created the mesh as shown in earlier posts in this series, we began work on building the 0 (time) folder.
The procedure was to copy all the existing files from the current running case, and then to examine each one to see if must be adjusted for the new mesh.
Here is the Co file in the 0 folder:
cat 0/Co
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format binary;
class volScalarField;
location "0";
object Co;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //dimensions [];
internalField uniform 0;
boundaryField
{
wallSection
{
type zeroGradient;
}
frontAndBackPlanes
{
type empty;
}
heatedWall
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type zeroGradient;
}
}// ************************************************************************* //
And here is ChatGPT5's evaluation of the file in the context of our new mesh:
Yes — those Co settings are perfectly consistent with your new mesh as long as the patch names match (and they do: wallSection, frontAndBackPlanes, heatedWall, inlet, outlet).
A couple of practical notes:
Co is a diagnostic field (dimensionless Courant number). It doesn’t drive the physics the way U, p, T, etc. do. So having:
internalField uniform 0;
zeroGradient on most patches
empty on frontAndBackPlanes
…is totally fine.
For the 2D-extruded mesh, frontAndBackPlanes { type empty; } is the key requirement, and you’ve got it.
In many cases, you don’t even need a 0/Co file (OpenFOAM can compute/report Courant without storing a Co volScalarField), but keeping it doesn’t hurt if your workflow expects it.
The next file we will examine is p (pressure)
(th)
Offline
Like button can go here
#243 Yesterday 10:11:47
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
Here is the p file in the 0 folder (p is for pressure - 0 is the default time folder)
cat 0/p
/*--------------------------------*- 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 volScalarField;
location "0";
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //dimensions [ 1 -1 -2 0 0 0 0 ];
internalField uniform 101325;
boundaryField
{
wallSection
{
type zeroGradient;
}
frontAndBackPlanes
{
type empty;
}
heatedWall
{
type zeroGradient;
}
inlet
{
type fixedFluxPressure;
}
outlet
{
type waveTransmissive;
field p;
gamma 1.4;
value uniform 20000;
}
}// ************************************************************************* //
And here is ChatGPT5's evaluation of p...
Yes — your p file is basically consistent with the new “simple pipe” mesh (same patch names, 2-D with frontAndBackPlanes as empty, etc.), and your two key modeling choices also make sense:
Preload at 1 bar (internalField uniform 101325;): totally reasonable. It’s a convenient “start filled” initial condition so the solver has something physical everywhere on step 1. You can keep it.
Outlet set to a low “vacuum-ish” back pressure (value uniform 20000;) with waveTransmissive: also reasonable, and it matches what you’ve been doing to avoid overly stiff numerics.
That said, there are a couple of important details to watch with this exact file:
1) inlet with fixedFluxPressure
This is not automatically “wrong,” but it’s only correct if your velocity/mass-flow boundary condition is imposed in U/phi in a compatible way.
fixedFluxPressure is typically used when the flux (phi) is being prescribed by the velocity boundary condition (or a mass flow model), and pressure is adjusted consistently to satisfy continuity.
If you later change inlet strategy (e.g., pressure-driven inlet), you might swap it — but for the “mass-flow inlet + compatible pressure” style you’ve been using, it’s consistent.
Rule of thumb: if your U file at the inlet is still something like a mass-flow / flowRate style condition (or otherwise constraining flux), then fixedFluxPressure is a good pairing.
2) outlet as waveTransmissive is good — but your gamma may be better taken from thermo
You set:
gamma 1.4; (you’ve used that as a workable approximation)
For hydrogen, gamma is not always 1.4 (especially as temperature rises), but using 1.4 as a stabilizing approximation is fine for now. If later you want to tighten physical correctness, we can look at whether your solver/thermo setup supports a better coupling.
3) The 2-D constraint is correctly applied
You have:
frontAndBackPlanes { type empty; }
That matches what checkMesh reported: Mesh has 2 solution directions (1 1 0).
Just make sure U, T, e (or h), rho, etc. also have frontAndBackPlanes set to empty where required (especially the vector fields).
4) About the “vacuum” value
Your outlet fallback:
value uniform 20000; (20 kPa)
This is a good “numerical vacuum” choice. If you make it too small (like 0 Pa), you can make the compressible system stiffer and sometimes trigger ugly behavior. So: keep 20 kPa unless you have a reason to sweep it.
What I’d check next (quick sanity list)
To confirm the p file will behave the way you intend in the new mesh:
0/U at inlet and outlet: does it still impose the inlet mass flow the way you want?
Thermo energy field: are you solving T + e, or h, etc.? (Make sure p BCs are compatible with that formulation.)
0/phi and controlDict massFlow function objects: confirm the sign convention is still as expected (you saw massFlow_in = -2, which is normal with OpenFOAM’s patch normal conventions).
(th)
Offline
Like button can go here
#244 Yesterday 10:14:17
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
In this post, ChatGPT5.2 recapitulates the two files we've seen so far, and looks at the big picture:
Yes — both the Co and p files you’ve shown are consistent with your new blockMesh-based pipe mesh, and they align well with the way you’ve been running extendedMerlin. I’ll walk through them carefully and point out the few things to keep in mind, without asking you to change anything prematurely.
1. Co field — Correct and appropriate
dimensions [];internalField uniform 0;
✔️ Correct.
You’re not actively transporting Courant number; it’s a diagnostic field written by function objects / solvers. Starting at zero is standard.Boundary conditions
wallSection zeroGradient;
heatedWall zeroGradient;
inlet zeroGradient;
outlet zeroGradient;
frontAndBackPlanes empty;✔️ All of these are appropriate and consistent with:
A 2-D extrusion mesh (your frontAndBackPlanes = empty)
A pipe-like domain
No artificial constraints on Co evolution
There is no conflict between this Co file and your new mesh.
✅ Verdict: leave Co exactly as-is.
2. p field — Also correct, with one important conceptual confirmation
internalField uniform 101325;✔️ This is exactly right given your stated intent:
You are preloading hydrogen at 1 bar
You want a clean zero-time startup
You want to avoid numerical shock from vacuum initialization
This is a deliberate modeling choice, and a good one for your goals.
Wall and symmetry boundaries
wallSection
{
type zeroGradient;
}heatedWall
{
type zeroGradient;
}frontAndBackPlanes
{
type empty;
}✔️ All correct.
Pressure at solid walls should not be prescribed
zeroGradient is the right choice
empty is mandatory for 2-D
Inlet
inlet
{
type fixedFluxPressure;
}✔️ Correct given your solver setup.
This pairs properly with:
Mass-flow-controlled inlet via U
Compressible formulation
Smooth pressure adjustment to enforce flux
You’ve already proven this works in production runs.
Outlet (this is the critical one)
outlet
{
type waveTransmissive;
field p;
gamma 1.4;
value uniform 20000;
}✔️ This is entirely consistent with what you’ve been doing and what you want now.
Key points:
waveTransmissive is the correct choice to minimize reflections
value uniform 20000 is a numerically sane “near-vacuum”
You are explicitly avoiding true vacuum stiffness — wisely
This outlet condition is mesh-agnostic, so it carries over cleanly to your new pipe mesh.
3. Big-picture confirmation (important)
What you’ve done here is actually exactly right for the transition you’re making:
Same patch names ✅
Same physics files ✅
Same BC philosophy ✅
New mesh geometry only ✅
That means:
If something behaves differently in the new mesh, it will be because of geometry — not configuration drift.
That’s the gold standard for a controlled CFD experiment.
4. One gentle forward-looking note (not a change request)
When you eventually:
Remove the washer
Add a real nozzle
Eliminate the enclosing “box”
You may later experiment with:
Lower outlet value (e.g. 5–10 kPa)
Or a slightly different non-reflecting outlet formulation
But now is not that time.
For validating:
Heating
Mass flow
Stability
Thrust accounting
? Your current p file is correct and defensible.
Final verdict
✔️ Yes — these settings are fully consistent with your new mesh.
✔️ You should proceed exactly as planned.
✔️ Documenting this as a tutorial is an excellent idea — this is a rare, clean example of “physics preserved, geometry swapped.”
(th)
Offline
Like button can go here
#245 Yesterday 10:18:56
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
The next file we will study is T (temperature)
cat 0/T
/*--------------------------------*- 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 volScalarField;
location "0";
object T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //dimensions [ 0 0 0 1 0 0 0 ];
internalField uniform 21; // was 300
boundaryField
{
wallSection
{
type zeroGradient;
}
frontAndBackPlanes
{
type empty;
}
heatedWall
{
type zeroGradient;
}
inlet
{
type fixedValue;
value uniform 21; // was 300
}
outlet
{
type zeroGradient;
}
}// ************************************************************************* //
Here is ChatGPT5.2's evaluation:
Yep — your 0/T file looks internally consistent with the new pipe mesh and with the way you’ve been running the case.
What looks good
internalField uniform 21;
Perfectly fine (and far more physically meaningful than the “training wheels 300 K”), as long as your thermo model supports that temperature range for H₂ (it should).frontAndBackPlanes = empty
Correct for your thin-extruded 2D mesh.inlet = fixedValue 21
This matches your intent: inject cold H₂ at the inlet.walls (wallSection, heatedWall) = zeroGradient
This is consistent with an adiabatic wall temperature boundary condition (i.e., you’re not prescribing wall temperature). Since your heating is volumetric via fvModels, you typically do want wall T to be free to float, so zeroGradient is reasonable.The one thing I’d consider revisiting (later)
outlet = zeroGradient for T
This is common and often OK, especially when the outlet is mostly outflow.
But if you ever see backflow at the outlet (which can happen during transients, pulses, or reflections), zeroGradient can allow unhelpful temperature values to “enter” from the outlet numerically.If you want a slightly more robust outlet for transients, you can use an “inletOutlet”-style temperature BC (name varies by OpenFOAM flavor/version), i.e.:
outflow: zeroGradient
inflow (reverse flow): fixedValue to something sane (often your inlet temperature or ambient)
If you tell me whether you’re on OpenFOAM Foundation v12 (it looks like it) and what inletOutlet variants you have available, I can give you the exact outlet stanza.
But as-is: your T file is a solid, sensible starting point for the new mesh.
(th)
Offline
Like button can go here
#246 Yesterday 11:01:33
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
In this post we will look at the U (velocity) file. In this one, we are prescribing the ramp for input.
cat 0/U
/*--------------------------------*- 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 volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //// Dropped feed to 1 kg/s 2026/01/06 for run to 16
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
wallSection
{
type fixedValue;
value uniform (0 0 0);
}
frontAndBackPlanes
{
type empty;
}
heatedWall
{
type fixedValue;
value uniform (0 0 0);
}
inlet
{
type flowRateInletVelocity;
massFlowRate table
( // seconds kilograms
(0.00 0.00) // kg/s
(0.50 0.20)
(1.00 1.00)
(2.00 2.00) // reach target more gently
(10000 2.00) // hold
);
rho rho;
value uniform (0 0 0);
}
outlet
{
type pressureInletOutletVelocity;
value uniform (0 0 0);
}
}
Here is ChatGPT5.2's assessment:
Please note that this reply contains an example of ChatGPT5.2 not being able to remember details from a long running project. In this case, it did not remember that we have been running simulations with three different flow rates: 1 kg/s, 2 kg/s and 20.
I deal with this limitation by refreshing data as necessary. In this case a refresh is not needed because the comment was advisory. I composed a reply reminding ChatGPT5.2 that we will be testing with multiple feed rates because we want to be able to throttle the engine depending upon navigation requirements.
Lower flow rate will increase ISP so we will be using propellant more efficiently. On the other hand, higher flow rate will allow greater thrust at the cost of efficiency.
Thought for 12s
Yep — structurally, that U file is compatible with your new mesh as long as the patch names match (they do: wallSection, frontAndBackPlanes, heatedWall, inlet, outlet) and your inlet/outlet patch normals are sensible.
That said, I see two important gotchas and one “vacuum” note:
1) Your inlet table is not actually 1 kg/s
You’ve got comments saying “dropped feed to 1 kg/s” and “reach target more gently”… but the table ramps to 2.00 kg/s at t=2.00 and holds it:
(1.00 1.00)
(2.00 2.00)
(10000 2.00)So: this is a 2 kg/s schedule, not 1 kg/s.
If you want the “gentle” ramp to 1 kg/s, a clean version is:
massFlowRate table
(
(0.00 0.00)
(0.50 0.20)
(1.00 1.00)
(10000 1.00)
);If you want 2 kg/s but still gentle, keep yours (it’s fine), just fix the comment so future-you doesn’t get tricked.
2) flowRateInletVelocity will push flow along the inlet patch normal
That’s usually exactly what you want: OpenFOAM computes a uniform (or patch-distributed) inlet velocity so that the mass flow rate is satisfied.
But it assumes the patch normal points into the domain. If your inlet patch normal is flipped, you’ll get the “right” magnitude but the wrong sign/direction weirdness.
Quick sanity checks:
In the log, you’ve been seeing sum("inlet") of phi = -2 for 2 kg/s. That’s consistent with a common sign convention (inlet negative). So don’t panic about the sign by itself.
If you want to verify direction, in ParaView you can show U vectors near the inlet and confirm they point down the pipe.
3) Outlet “vacuum”: why we keep a small value
You’re remembering correctly: true deep vacuum values are numerically brutal in compressible solvers with wave-type outlet conditions, especially with any chance of transient reverse flow.
LEO “pressure” is effectively near-zero for CFD purposes (orders of magnitude below 1 Pa), but using ~0 Pa as a boundary target tends to make the outlet too stiff and can trigger instability.
Your current p outlet fallback value uniform 20000; is 20 kPa, which is not vacuum physically — it’s a numerical stabilizer.
If you want to push it “more vacuum-like” carefully, a good approach is to step it down across runs:
try value uniform 5000; (5 kPa),
then 2000;,
then 1000;,
and only later consider 100; if it stays stable.
That gets you closer to “vacuum behavior” without instantly destabilizing the run.
Bottom line on your U file
Walls: fixedValue (0 0 0) on wallSection and heatedWall is correct (no-slip).
frontAndBackPlanes: empty is correct for your thin-Z 2D extrusion.
inlet: correct if the table matches your intended mass flow (fix that 1 vs 2 mismatch).
outlet: pressureInletOutletVelocity is reasonable and pairs fine with your pressure outlet strategy.
If you paste your new mesh’s blockMeshDict patch definitions (just the patch blocks), I can tell you immediately whether the inlet/outlet normals and axis assumptions line up with your intended flow direction.
The next file we will examine is probably not needed for the current case. It may be left over from a time when we were using convective heating.
Our current heating mode is called 'volume" heating. It mathematically simulates perfect operation of our 240 meter heat pipe. 40 MW of thermal energy is mathematically distributed over 24,000 cells in our heating subsystem.
(th)
Offline
Like button can go here
#247 Yesterday 12:27:59
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
This post is about a file I copied over from the running case. This file may not be required in the current case, because we are now using "volume" heating instead of convection heating.
cat 0/wallHeatFlux
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 12
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format binary;
class volScalarField;
location "0";
object wallHeatFlux;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //dimensions [1 0 -3 0 0 0 0];
internalField uniform 0;
boundaryField
{
wallSection
{
type calculated;
value uniform 0;
}
frontAndBackPlanes
{
type empty;
}
heatedWall
{
type calculated;
value uniform 0;
}
inlet
{
type calculated;
value uniform 0;
}
outlet
{
type calculated;
value uniform 0;
}
}
As I hoped, ChatGPT5.2 recognized this file as not needed for our new mesh. I deleted the file.
The 0 directory is now ready for a run. However, we still have the system and constant folders to examine.
Here is the assessment for the unneeded file:
Yes — your 0/wallHeatFlux file is almost certainly not needed for the new mesh + volume heating (fvModels) path, and in many cases it’s best to remove it to avoid confusion.
Why it exists
wallHeatFlux is typically a derived/calculated field written by some solvers/functionObjects when you’re doing wall-based heat transfer (e.g., a wall heat flux boundary condition, or function objects that compute wall heat flux). Your file even says:
type calculated;
internalField uniform 0;
That strongly suggests it’s not an input you “set”—it’s an output/diagnostic.
What I recommend for your new “pipe-only” mesh case
If you are only using volumetric heating (your constant/fvModels heaterCells source) and not using wall heating BCs or a functionObject that writes wallHeatFlux:
✅ Delete 0/wallHeatFlux from the new case.
If the solver/functionObject wants it later, it will typically recreate it.
When you would keep it
Keep (or allow it to be generated) only if you are doing either of these:
You have wall heat transfer physics enabled that needs wall heat flux for post-processing.
You have a functionObject that explicitly computes/writes it and you want that diagnostic.
Quick sanity check (fast)
In your new case, search system/controlDict for anything like:
wallHeatFlux
heatFlux
wallHeat
wallHeatTransfer
If nothing references it, it’s safe to remove.
If you paste your system/controlDict functionObjects section (or the whole file), I can tell you definitively whether anything in your setup will try to write/use wallHeatFlux.
(th)
Offline
Like button can go here
#248 Yesterday 14:36:22
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
We'll take a break here, from working on the new mesh, to take a look at the latest set of sprints for extendedMerlin.
Here is the end of the tail from the run after reaching 8 seconds:
ExecutionTime = 98218.89 s ClockTime = 98540 s
volFieldValue p_avg write:
volAverage(heaterCells) of p = 334953.3volFieldValue T_avg write:
volAverage(heaterCells) of T = 716.2811volFieldValue rho_mass write:
volIntegrate(heaterCells) of rho = 3.011971volFieldValue T_min_global write:
min() of T = 22.45688 at location (-240.275 0.16 0) in cell 37871volFieldValue T_max_global write:
max() of T = 1543.261 at location (-0.182823 -0.1878 3.407286e-18) in cell 3401surfaceFieldValue massFlow_in write:
sum("inlet") of phi = -2surfaceFieldValue massFlow_out write:
sum("outlet") of phi = 1.979579End
Follow up plans include:
1) Thrust plot
2) Plots generated from the logs over all 8 seconds
3) Any comments ChatGPT5.2 may have
4) Any comments NewMars members may have.
As a reminder, we ran this case with the washer installed in the throat. We fed 2 kg/s and heated with 40 MW.
The purpose of the run beyond 4 seconds was to see if the system settled. It was still changing at 4 seconds.
No surprises. Fluctuations continued due to the flat washer in the intake.
The ** real ** information will come from the 8 plots to be generated from the logs.
The 8 plots are 7 ... on plot had no data. And not much dramatic showed up. The system has not settled but changes are minor in the later seconds.
Pressure holds steady above 3 bar
Mass is still accumulating in the heater pipe but almost imperceptibly
Average temperature is leveling off
Maximum Temperature is surging with fluctuations seen in the thrust plot
Minimum Temperature is holding steady at the intake as expected.
Time slices are spiking with temperature due to fluctuations
The control for time for equations reaches the prescribed maximum and stays there.
It seems to me we've wrung about as much value as we can out of this test system. The "washer" was useful for getting a hint about ISP but it is time to move on to the new mesh. The new mesh will be installed in this case when it is working in foamWorkroom.
(th)
Offline
Like button can go here
#249 Yesterday 15:38:07
- SpaceNut
- Administrator
- From: New Hampshire
- Registered: 2004-07-22
- Posts: 30,224
Re: OpenFOAM
After thinking about the thermal engine type and the last google meeting as GW put it the oscillation results from the modeling.
similar to the engine is ION and NTP which creates thrust from temperature change.
Ion engine thrust oscillation, particularly regarding instability, often manifests as low-frequency "breathing modes" (typically 10–30 kHz) in Hall thrusters due to oscillating plasma and neutral densities.
Based on modeling and experimental data:Initial Transient Period: Thrust instability is most pronounced shortly after ignition as the thruster reaches thermal equilibrium, with significant noise observed between 150–300 seconds.
3-Second Behavior: While the major thermal stabilization occurs over hundreds of seconds, high-voltage breakdown in grid-based ion thrusters can cause immediate, short-duration thrust fluctuations, potentially lasting only a few seconds.
Modeling Approach: A zero-dimensional (0D) "predator-prey" model using two coupled ordinary differential equations is commonly used to characterize these oscillations in plasma and neutral density.
Stabilization Mechanisms: The oscillation behavior is highly dependent on the stability of the beam neutralization. If the neutralizer fails to provide electrons to neutralize the ion beam, the thruster becomes unstable, causing significant, rapid oscillations in thrust.
For modern, well-designed ion engines, once the initial startup sequence is complete (well after the first 3 seconds), the thrust is generally stable, operating with a very low-level, high-frequency noise
Thrust oscillation in hydrogen-based thermal rocket engines, particularly Nuclear Thermal Propulsion (NTP) or hydrogen-fed thrusters, is a critical transient behavior that often initiates during startup, typically appearing within the first few seconds of operation (3–40 seconds) as the engine approaches steady-state chamber pressure.
Key Aspects of the Model:
Timeframe: Instabilities often manifest immediately following ignition (within 3 seconds) or during the initial ramp-up of hydrogen flow.Mechanism: The oscillations are often driven by acoustic resonance in the combustion chamber (or reactor core) coupling with the propellant feed system, particularly when dealing with cryogenic hydrogen.
Oscillation Characteristics: These are often low-frequency (0–300 Hz) instabilities, characterized by periodic pressure, heat, and flow rate fluctuations that can lead to significant thrust oscillations.
Modeling and Simulation: Approach: Models typically use a transient analysis (e.g., using FEM or CFD in software like ANSYS) to simulate the first 3 seconds of ignition and subsequent startup phase.
Instability Drivers: In \(H_{2}/O_{2}\) systems, low-temperature hydrogen (e.g., 95 K) is known to induce self-excited instabilities. In NTP systems, the rapid heating of hydrogen to over 2500 K through the core can cause pressure spikes.
Frequency Range: Simulations frequently focus on the first longitudinal resonance (1-L) of the injector/combustion chamber. Impact at 3
Seconds:After 3 seconds, the engine is usually still in its ignition transient phase. High-frequency, self-sustained oscillations can create pressure waves with amplitudes of 80% or more of the mean chamber pressure, leading to significant thermal, pressure, and structural loads on the chamber, particularly at the throat.
Offline
Like button can go here
#250 Yesterday 19:13:41
- tahanson43206
- Moderator
- Registered: 2018-04-27
- Posts: 23,880
Re: OpenFOAM
With this post we resume our tutorial on development of the new mesh for the Optical Plane vessel engine.
As we resume, we have set up the 0 (Zero) folder and the polyMesh folder in the constant folder.
ChatGPT5.2 looked at the file listing for the constant folder and identified only fvModels as needing attention. This is the file that controls operation of the volume heating feature that allows us to feed 40 MW of energy into the heat pipe. We need to make sure the zone to be heated in the new mesh matches the zone that we heated in extendedMerlin. I am hoping no changes are needed, but the face numbers and point numbers likely changed, so we may have to make adjustments.
Here is the list of files in the constant folder after they were copied from the running case:
ls constant
fvModels physicalProperties polyMesh thermophysicalProperties turbulenceProperties
polyMesh is a folder containing the mesh created as shown in posts above this one.
Here is ChatGPT5.2's evaluation:
Great — that constant/ inventory is exactly what we want at this stage, and yes: keeping the same patch names (inlet outlet wallSection heatedWall frontAndBackPlanes) is what protects all the work you’ve already done in 0/, system/, post-processing, thrustCalc, etc.
What to keep in constant/ for the new mesh case
1) constant/polyMesh/Keep it (obviously) — this is the new mesh.
Two quick sanity checks that matter a lot:
constant/polyMesh/boundary patch list must include exactly your 5 patches (names spelled identically).
Patch types should match the intent:
frontAndBackPlanes → empty (since you’re doing a thin quasi-2D model)
wallSection and heatedWall → typically wall
inlet, outlet → patch
If those match, you’ve preserved compatibility with nearly everything else.
2) constant/thermophysicalProperties
Keep it as-is from the working case. This is where your hydrogen model (perfect gas vs real gas, Cp model, transport, etc.) lives. A mismatch here will absolutely destabilize a case even if everything else is perfect.
3) constant/turbulenceProperties
Keep it. Even if you’re running laminar, this file is still part of the “known-good” configuration.
4) constant/physicalProperties
Keep it. Depending on your solver and setup, this may be required (or harmless). Since it came from the working run, don’t touch it yet.
5) constant/fvModels
Keep it and this is the big win: having fvModels in constant/ is consistent with the working volume-heating approach you converged on.
One caution: the contents of fvModels usually reference:
a region selection (cellZone / cellSet / topoSet results), and/or
a list of cells (like heaterCells)
Since your new mesh is different, the zone/cellSet names can stay the same, but you’ll need to recreate the selection for the new mesh (so the heating applies to the correct cells). That’s a later step, but it’s the one place where “same patch names” alone isn’t enough.
Next best file to review
If you paste:
constant/fvModels
and (later) whatever you’re using to define the heater region (often system/topoSetDict or a topoSetDict.* plus your cellSet/cellZone workflow)
…I can tell you exactly what must be re-targeted for the new mesh, and what can remain unchanged.
We'll work on fvModels in the next post.
(th)
Offline
Like button can go here