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,904

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,904

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,904

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 Yesterday 07:58:05

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

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 Yesterday 08:12:40

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

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 Yesterday 11:17:23

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

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 Yesterday 13:30:34

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

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 Yesterday 18:50:38

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

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 Yesterday 20:28:49

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

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 Yesterday 21:16:38

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

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 Yesterday 22:04:31

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

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 Yesterday 22:12:46

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

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

Board footer

Powered by FluxBB