Debug: Database connection successful OpenFOAM (Page 11) / Science, Technology, and Astronomy / New Mars Forums

New Mars Forums

Official discussion forum of The Mars Society plus New Mars Image Server

You are not logged in.

Announcement

Announcement: This forum is accepting new registrations via email. Please see Recruiting Topic for additional information. Write newmarsmember[at_symbol]gmail.com.

#251 2026-01-17 19:15:10

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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

#252 2026-01-17 22:00:34

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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

#253 2026-01-17 22:07:25

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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.4

volFieldValue T_avg write:
    volAverage(heaterCells) of T = 311.5977

volFieldValue rho_mass write:
    volIntegrate(heaterCells) of rho = 3.829058

volFieldValue T_min_global write:
    min() of T = 1.194984 at location (4.249592 -1.607025 4.607311e-18) in cell 13793

volFieldValue T_max_global write:
    max() of T = 338.8237 at location (-0.182823 -0.1878 3.407286e-18) in cell 3401

surfaceFieldValue massFlow_in write:
    sum("inlet") of phi = -1.000001

surfaceFieldValue massFlow_out write:
    sum("outlet") of phi = 9.772167

End

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

#254 2026-01-18 07:58:05

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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

#255 2026-01-18 08:12:40

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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

#256 2026-01-18 11:17:23

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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

#257 2026-01-18 13:30:34

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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

#258 2026-01-18 18:50:38

SpaceNut
Administrator
From: New Hampshire
Registered: 2004-07-22
Posts: 30,284

Re: OpenFOAM

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 / Δh

 Step-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_e

 BBcode 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

#259 2026-01-18 20:28:49

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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 time

Create 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 present

Courant 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 mode



Starting 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.8


surfaceFieldValue massFlow_in: <<
    total faces  = 40
    total area   = 0.04


surfaceFieldValue massFlow_out: <<
    total faces  = 40
    total area   = 0.04


volFieldValue 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 = 0

surfaceFieldValue massFlow_out write:
    sum("outlet") of phi = 0

Courant Number mean: 0 max: 0
deltaT = 1.2e-06
Time = 1.2e-06s

diagonal:  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

#260 2026-01-18 21:16:38

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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 heating

volFieldValue T_avg write:
    volAverage(heaterCells) of T = 402.8296 << We are definitely seeing heating because we started at 21 Kelvin

volFieldValue rho_mass write:
    volIntegrate(heaterCells) of rho = 3.048076 << Mass is presumably falling from the initial load

volFieldValue 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 expected

volFieldValue 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 exit

surfaceFieldValue 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 seconds

surfaceFieldValue massFlow_out write:
    sum("outlet") of phi = 4.657116 << We are seeing a high value because we are flushing the initial load

End << 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

#261 2026-01-18 22:04:31

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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-correct

From 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

#262 2026-01-18 22:12:46

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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 bar

volFieldValue T_avg write:
    volAverage(heaterCells) of T = 1459.301 << Heating is happening in a smaller volume 9600 vs 24000

volFieldValue rho_mass write:
    volIntegrate(heaterCells) of rho = 2.155162 << flush is proceeding nicely while we ramp up to 2 kg/s

volFieldValue 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 Kelvin

volFieldValue 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 exit

Courant 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

#263 2026-01-19 08:09:45

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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 over

volFieldValue T_avg write:
    volAverage(heaterCells) of T = 897.7812 << There is substantial heating going on

volFieldValue rho_mass write:
    volIntegrate(heaterCells) of rho = 4.952744 << That mass figure seems high to me. We have 2 kg/s coming in

volFieldValue 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 Kelvin

volFieldValue 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 Kelvin

surfaceFieldValue massFlow_in write:
    sum("inlet") of phi = -1.999996 << We are feeding 2 kg/s into the intake

surfaceFieldValue massFlow_out write:
    sum("outlet") of phi = 1.773398 << Output should eventually match input. We are not there yet

End

We'll pause here to create plots, and run ParaView to see the animation of the run.

(th)

Offline

Like button can go here

#264 2026-01-19 14:19:40

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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.

file.php?id=204

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

#265 2026-01-19 18:41:30

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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 pressure

  • Ambient 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 tons

  • Thrust 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 momentum

  • Reaction 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 nozzle

Pressure 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/s

Interpretation:

  • []System has not yet purged all preload gas
    []Some of the original 1-bar hydrogen inventory is still exiting

  • Pressure 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 kN

  • Then slowly declines toward ~29 kN by 4 s

This suggests:

  • []Initial acceleration phase as pressure builds
    []Partial mass balance stabilization

  • Still 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 cm

  • Planned throat area ≈ 0.0007 m²

That is roughly a 60× reduction in throat area.

This will:

  • []Raise chamber pressure substantially
    []Drive the flow choked

  • Greatly 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 input

The 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 throat

  • Define 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

#266 2026-01-19 18:48:37

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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.9

volFieldValue T_avg write:
    volAverage(heaterCells) of T = 757.2487

volFieldValue rho_mass write:
    volIntegrate(heaterCells) of rho = 4.465979

volFieldValue T_min_global write:
    min() of T = 26.83986 at location (-239.5 0.175 2.292485e-17) in cell 8880

volFieldValue T_max_global write:
    max() of T = 1789.563 at location (-0.5 -0.195 -2.162287e-17) in cell 239

surfaceFieldValue massFlow_in write:
    sum("inlet") of phi = -2

surfaceFieldValue massFlow_out write:
    sum("outlet") of phi = 1.893913

End

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

#267 2026-01-20 08:01:24

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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 9s

You’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

#268 2026-01-20 08:39:58

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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 bar

volFieldValue T_avg write:
    volAverage(heaterCells) of T = 757.2487 << This is average over 9600 cells

volFieldValue 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 door

volFieldValue 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 3000

surfaceFieldValue massFlow_in write: << We are pushing 2 kg/s in at the intake
    sum("inlet") of phi = -2

surfaceFieldValue massFlow_out write: << and we are NOT YET at equilibrium after 16 seconds
    sum("outlet") of phi = 1.893913

End

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 result is literally ** unbelievable ** so please do NOT believe it.  I don't know what went wrong but I suspect the thrust calculator did not make the transition from extendedMerlin to foamWorkroom intact.  With that caution in mind, here is the plot after 16 seconds:

file.php?id=205

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

#269 2026-01-20 10:47:21

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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.

file.php?id=205
Thrust (unbelievable)

file.php?id=206
Pressure - does this make sense?

file.php?id=207
Mass - looks reasonable given 2 kg/s input

file.php?id=208
Temperature average over 9600 cells

file.php?id=210
Minimum temperature near intake - looks reasonable

file.php?id=209
Maximum temperature - note 5000 max cut off by safety settings

file.php?id=212
The time slices held study ... the equations were loafing.

file.php?id=211
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

#270 2026-01-20 13:11:00

GW Johnson
Member
From: McGregor, Texas USA
Registered: 2011-12-04
Posts: 6,126
Website

Re: OpenFOAM

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

#271 2026-01-20 18:46:40

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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 exitVelocity

This is correct only if:

  • [] The outlet is a real exhaust plane
    [] Momentum is carried away downstream

  • The 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 vacuum

In 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 surfaces

  • The 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 patch

So it assumes:
"This must be thrust."

But in your case:

  • [] The outlet plane is not a nozzle exit
    [] The momentum is not free-stream exhaust

  • The 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 N

The 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 massFlowIn

This is exactly what your logs show:

  • [] inlet ~2.0 kg/s
    [] outlet ~1.77 kg/s

  • pressure 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 collimated

  • No 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 velocity

  • Redirects lateral pressure into axial momentum

Only then does:

  • [] Pressure drop downstream
    [] Exhaust carry momentum away

  • thrustCalc 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 nozzle

  • The 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 nozzle

  • Define 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

#272 2026-01-20 19:22:22

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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 patch

In 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 area

  • Mass flow imposed at the inlet

That means:

  • []There is no proper jet formation
    []There is no supersonic expansion

  • There 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 boundary

  • And that force being interpreted as vehicle thrust

This is physically equivalent to:

  • []A pressure-fed piston
    []Not a rocket exhaust

Why 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 thrust

This is allowed — but it is not propulsion.

A rocket requires:

  • []Momentum leaving the system
    []At high directed velocity

Your 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 increases

  • Mass 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 work

  • That 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 irrelevant

In your case:

  • []Exit velocity is modest
    []Exit momentum is too small for 28 kN

  • Pressure 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 consistent

But it also proves:

  • []A nozzle (or expansion surface) is absolutely required
    []Without it, thrust numbers are not meaningful

What changes once a nozzle is added

When you add:

  • []A throat
    []An expansion wedge or bell

Then:

  • []Pressure converts to directed velocity
    []Momentum leaves the system
    []Inlet force becomes negligible
    []Thrust becomes physically real

At 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 momentum

You’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

#273 2026-01-20 19:55:52

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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

#274 2026-01-20 21:47:49

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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 11s

Yes — 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

#275 Yesterday 05:37:38

tahanson43206
Moderator
Registered: 2018-04-27
Posts: 23,959

Re: OpenFOAM

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 increasing

volFieldValue T_avg write:
    volAverage(heaterCells) of T = 729.9388

volFieldValue rho_mass write:
    volIntegrate(heaterCells) of rho = 5.327563 << mass still increasing

volFieldValue T_min_global write:
    min() of T = 26.83269 at location (-239.5 0.175 2.292485e-17) in cell 8880 << Temperature at intake steady

volFieldValue T_max_global write:
    max() of T = 1703.316 at location (-0.5 -0.195 -2.162287e-17) in cell 239 << Temperature at exit

surfaceFieldValue massFlow_in write: << mass flow in steady
    sum("inlet") of phi = -2

surfaceFieldValue massFlow_out write:
    sum("outlet") of phi = 1.972983 << Still below intake after 32 seconds but close

End

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:
file.php?id=214

(th)

Offline

Like button can go here

Board footer

Powered by FluxBB