OpenFOAM 14, caso 01: flusso laminare tra lastre piane validato con Poiseuille

Questo è il primo caso del percorso Imparare OpenFOAM 14: un flusso laminare, stazionario e incomprimibile tra due lastre piane. L’ho scelto apposta perché ha una soluzione esatta, il profilo di Poiseuille, e quindi il risultato numerico si può controllare con i numeri e non solo a occhio. Alla fine il profilo calcolato coincide con quello teorico entro lo 0,2%.

Quella che segue è la procedura completa e riproducibile del caso, con la validazione analitica e la visualizzazione in ParaView, nello stesso ordine della dispensa originale: OpenFOAM 14 (Foundation / CFD Direct) e ParaView 5.13.2.

1. Scopo e risultato atteso

Questa dispensa descrive la procedura completa necessaria per costruire, eseguire, validare e visualizzare un flusso laminare incomprimibile, stazionario e bidimensionale tra due lastre piane parallele con condizione di aderenza (no-slip), usando OpenFOAM 14.

Il caso è volutamente semplice perché ammette una soluzione di Poiseuille in forma chiusa. Il risultato numerico può quindi essere verificato in modo quantitativo e non solo visivo.

Principali risultati attesi:

Grandezza	Valore
Versione di OpenFOAM	14
Modulo solutore	incompressibleFluid tramite foamRun
Mesh	200 × 40 × 1 = 8.000 celle esaedriche
Velocità media in ingresso	0,010 m/s
Viscosità cinematica	1,0·10⁻⁶ m²/s
Numero di Reynolds sul diametro idraulico	Re = 200
Etichetta dell’ultima iterazione stazionaria	2000
Ux massima campionata (CFD)	0,014962547 m/s
Ux massima analitica	0,015000000 m/s
RMSE normalizzato	0,150343%
Errore assoluto massimo / Umax	0,187193%
Visualizzazione	ParaView 5.13.2 in client-server con pvserver

La soluzione numerica accettata riproduce il profilo analitico di Poiseuille con un errore nettamente inferiore allo 0,2% sull’intera sezione campionata.

2. Prerequisiti

Devono essere disponibili i seguenti comandi:

foamVersion
command -v blockMesh
command -v checkMesh
command -v foamRun
command -v foamPostProcess
command -v pvserver

Versione di OpenFOAM attesa:

14

L’ambiente deve definire $FOAM_RUN. Il caso descritto qui viene creato in:

$FOAM_RUN/01_laminar_channel

3. Modello fisico

Dimensioni del canale:

Parametro	Valore
Lunghezza del canale L	0,500 m
Distanza tra le lastre H	0,010 m
Profondità numerica	0,001 m
Condizione davanti/dietro	empty

Fluido e condizioni di ingresso:

Parametro	Valore
Viscosità cinematica ν	1,0·10⁻⁶ m²/s
Velocità media in ingresso Umean	0,010 m/s

Il diametro idraulico per due lastre parallele molto larghe è circa

D_h = 2H = 0,020 m

e quindi

Re = U_mean·D_h / ν = (0,01 · 0,02) / 10⁻⁶ = 200.

Il flusso è volutamente ben all’interno del regime laminare.

4. Creare la cartella del caso

mkdir -p $FOAM_RUN/01_laminar_channel/{0,constant,system}
cd $FOAM_RUN/01_laminar_channel

5. Definizione della mesh

Creare system/blockMeshDict:

FoamFile
{
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}

scale 1;

vertices
(
    (0.0  -0.005  0.0)
    (0.5  -0.005  0.0)
    (0.5   0.005  0.0)
    (0.0   0.005  0.0)

    (0.0  -0.005  0.001)
    (0.5  -0.005  0.001)
    (0.5   0.005  0.001)
    (0.0   0.005  0.001)
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (200 40 1) simpleGrading (1 1 1)
);

edges
(
);

boundary
(
    inlet
    {
        type patch;
        faces
        (
            (0 4 7 3)
        );
    }

    outlet
    {
        type patch;
        faces
        (
            (1 2 6 5)
        );
    }

    walls
    {
        type wall;
        faces
        (
            (0 1 5 4)
            (3 7 6 2)
        );
    }

    frontAndBack
    {
        type empty;
        faces
        (
            (0 3 2 1)
            (4 5 6 7)
        );
    }
);

mergePatchPairs
(
);

La mesh strutturata contiene

200 × 40 × 1 = 8000

celle esaedriche.

6. Controllo del calcolo

Creare system/controlDict prima di lanciare le utility di mesh:

FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}

solver          incompressibleFluid;

startFrom       startTime;
startTime       0;

stopAt          endTime;
endTime         2000;

deltaT          1;

writeControl    runTime;
writeInterval   200;

purgeWrite      0;

writeFormat     ascii;
writePrecision  8;
writeCompression off;

timeFormat      general;
timePrecision   6;

runTimeModifiable true;

Poiché il caso usa steadyState in fvSchemes, i valori 1…2000 sono etichette di iterazione / pseudo-tempo, non un transitorio fisico di 2000 secondi.

7. Generare e verificare la mesh

blockMesh | tee log.blockMesh
checkMesh | tee log.checkMesh

Controllo utile:

grep -E "points:|faces:|cells:|hexahedra:|Mesh OK" log.checkMesh

Risultato atteso:

cells:          8000
hexahedra:      8000
Mesh OK.

8. Campi iniziali e condizioni al contorno

8.1 Campo di velocità: 0/U

FoamFile
{
    format      ascii;
    class       volVectorField;
    location    "0";
    object      U;
}

dimensions      [velocity];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.01 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    walls
    {
        type            noSlip;
    }

    frontAndBack
    {
        type            empty;
    }
}

8.2 Campo di pressione: 0/p

FoamFile
{
    format      ascii;
    class       volScalarField;
    location    "0";
    object      p;
}

dimensions      [kinematicPressure];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    walls
    {
        type            zeroGradient;
    }

    frontAndBack
    {
        type            empty;
    }
}

9. Proprietà del fluido e modello laminare

9.1 constant/physicalProperties

FoamFile
{
    format      ascii;
    class       dictionary;
    location    "constant";
    object      physicalProperties;
}

viscosityModel  constant;

nu              1e-06;

9.2 constant/momentumTransport

FoamFile
{
    format      ascii;
    class       dictionary;
    location    "constant";
    object      momentumTransport;
}

simulationType  laminar;

10. Schemi numerici

Creare system/fvSchemes:

FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}

ddtSchemes
{
    default         steadyState;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      bounded Gauss linearUpwind grad(U);
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

11. Solutori lineari e accoppiamento pressione-velocità

Creare system/fvSolution:

FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-10;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-10;
        relTol          0.05;
    }

    UFinal
    {
        $U;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor        no;
    nOuterCorrectors         1;
    nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
    fields
    {
        p   0.3;
    }

    equations
    {
        U   0.7;
    }
}

12. Eseguire il calcolo CFD

Eseguire un ultimo controllo della mesh e avviare il solutore:

checkMesh
foamRun | tee log.foamRun

Esaminare le ultime iterazioni:

tail -40 log.foamRun

Il calcolo validato ha raggiunto l’etichetta di iterazione finale 2000 e si è concluso normalmente con End.

Comportamento finale tipico:

Grandezza	Comportamento osservato
Tempo di esecuzione	circa 27,7 s
Residuo iniziale della pressione	ordine di 1e-10
Iterazioni della pressione	spesso 0–1, occasionalmente qualcuna in più
Errore di continuità locale	ordine di 1e-12
Terminazione	End normale

13. Soluzione analitica di Poiseuille

Per un flusso laminare completamente sviluppato tra lastre parallele, il profilo di velocità analitico è

u(y) = (3/2)·U_mean·[1 − (2y/H)²].

Con U_mean = 0,010 m/s e H = 0,010 m:

U_max,teoria = 1,5·U_mean = 0,015 m/s.

14. Post-processing del profilo di velocità

Creare system/graphUniform:

FoamFile
{
    format      ascii;
    class       dictionary;
    object      graphUniform;
}

start           (0.45 -0.004875 0.0005);
end             (0.45  0.004875 0.0005);
nPoints         40;

fields          (U);

axis            y;

#includeEtc "caseDicts/postProcessing/graphs/graphUniform.cfg"

Eseguire esplicitamente la function object:

foamPostProcess -solver incompressibleFluid -func graphUniform -latestTime

Individuare il risultato:

find postProcessing/graphUniform -type f

File atteso:

postProcessing/graphUniform/2000/line.xy

Esaminarlo:

head -10 postProcessing/graphUniform/2000/line.xy
tail -10 postProcessing/graphUniform/2000/line.xy
wc -l postProcessing/graphUniform/2000/line.xy

Il file contiene una riga di intestazione più 40 punti campionati. Le colonne sono:

y   U_x   U_y    U_z

Valori rappresentativi:

y [m]	Ux CFD [m/s]	Uy [m/s]
-0.004875	0.00073033707	-5.38e-17
-0.004125	0.0047752809	-2.79e-16
-0.003125	0.0091198502	-4.20e-16
-0.002875	0.010018726	-5.09e-16
+0.002875	0.010018726	-2.82e-15
+0.003875	0.0059737828	-4.36e-15
+0.004875	0.00073033707	-7.37e-16

La velocità trasversale è di fatto nulla e U_z è zero, come ci si aspetta in questo canale 2D.

15. Script di validazione quantitativa

Il comando seguente confronta ogni punto campionato con il profilo analitico di Poiseuille:

awk '
BEGIN {
    H=0.01
    Umean=0.01
    UmaxTheory=1.5*Umean
    umax=-1
}
$1 !~ /^#/ && NF >= 4 {
    y=$1
    u=$2

    uth=UmaxTheory*(1-(2*y/H)^2)

    e=u-uth
    ae=(e<0 ? -e : e)

    sumsq += e*e
    n++

    if (ae > maxerr) {
        maxerr=ae
        maxerr_y=y
        maxerr_cfd=u
        maxerr_theory=uth
    }

    if (u > umax) {
        umax=u
        umax_y=y
    }
}
END {
    rmse=sqrt(sumsq/n)

    printf "Points              : %d\n", n
    printf "CFD maximum Ux      : %.10f m/s\n", umax
    printf "Location of max     : %.7f m\n", umax_y
    printf "Theoretical Umax    : %.10f m/s\n", UmaxTheory
    printf "RMSE                : %.10e m/s\n", rmse
    printf "Normalized RMSE     : %.6f %%\n", 100*rmse/UmaxTheory
    printf "Maximum abs error   : %.10e m/s\n", maxerr
    printf "Max error / Umax    : %.6f %%\n", 100*maxerr/UmaxTheory
    printf "Max-error y         : %.7f m\n", maxerr_y
    printf "CFD there           : %.10f m/s\n", maxerr_cfd
    printf "Theory there        : %.10f m/s\n", maxerr_theory
}' postProcessing/graphUniform/2000/line.xy

Risultato validato:

Metrica	Valore
Punti confrontati	40
Ux massima CFD	0,0149625470 m/s
Posizione del massimo CFD	y = −0,0001250 m
Umax teorica	0,0150000000 m/s
RMSE	2,2551518927e-05 m/s
RMSE normalizzato	0,150343%
Errore assoluto massimo	2,8079000000e-05 m/s
Errore massimo / Umax	0,187193%
Posizione dell’errore massimo	y = +0,0001250 m
Valore CFD nel punto di errore massimo	0,0149625460 m/s
Teoria nello stesso punto	0,0149906250 m/s

Il campionamento usa 40 punti, quindi la mezzeria y = 0 non viene campionata esattamente. I campioni più vicini si trovano a circa y = ±0,000125 m.

16. Verifica della pressione

Per il flusso di Poiseuille piano, il gradiente di pressione cinematica atteso è

dp/dx = −12·ν·U_mean / H² = −1,2 × 10⁻³ m/s².

Su 0,5 m la corrispondente caduta di pressione cinematica è circa

6 × 10⁻⁴ m²/s²,

ovvero circa 0,6 Pa per l’acqua con ρ ≈ 1000 kg/m³.

Per campionare la pressione lungo la zona sviluppata, usare temporaneamente il seguente system/graphUniform:

FoamFile
{
    format      ascii;
    class       dictionary;
    object      graphUniform;
}

start           (0.30 0 0.0005);
end             (0.49 0 0.0005);
nPoints         40;

fields          (p);

axis            x;

#includeEtc "caseDicts/postProcessing/graphs/graphUniform.cfg"

Quindi eseguire:

foamPostProcess -solver incompressibleFluid -func graphUniform -latestTime

La pressione campionata deve diminuire in modo approssimativamente lineare nella zona completamente sviluppata.

17. Visualizzazione con ParaView in client-server

Per la visualizzazione interattiva si usa ParaView in modalità client-server.

17.1 Avviare il server

Controllare la versione installata:

command -v pvserver
pvserver --version

Versione validata:

paraview version 5.13.2

Creare il file marcatore di OpenFOAM:

cd $FOAM_RUN/01_laminar_channel
touch 01_laminar_channel.foam

Avviare pvserver:

pvserver --bind-address=0.0.0.0 --server-port=11111

17.2 Connettersi da ParaView sulla workstation

Usare la stessa versione di ParaView del server: 5.13.2.

In ParaView:

File -> Connect... -> Add Server

Configurare:

Name: CFD Server
Server Type: Client / Server
Host: <indirizzo IP del server>
Port: 11111

Dopo la connessione, File → Open naviga nel filesystem Linux remoto. Aprire:

$FOAM_RUN/01_laminar_channel/01_laminar_channel.foam

Selezionare internalMesh, abilitare U e p, quindi premere Apply.

17.3 Controlli visivi

All’etichetta di iterazione finale 2000:

1. Coloring → U → Magnitude
2. Rescale to Data Range
3. verificare una distribuzione di velocità simmetrica, con la velocità massima al centro del canale e velocità nulla alle pareti no-slip
4. Coloring → p
5. verificare una diminuzione regolare della pressione dall’ingresso all’uscita

Per una vista pulita delle linee di corrente, il canale può essere visualizzato come Outline oppure come superficie a bassa opacità.

18. Linee di corrente

Selezionare la sorgente OpenFOAM e creare:

Filters -> Stream Tracer

Impostazioni consigliate:

Proprietà	Valore
Vectors	U
Integration Direction	FORWARD
Seed Type	Line
Point1	(0.001, -0.0045, 0.0005)
Point2	(0.001, +0.0045, 0.0005)
Resolution	20
Maximum Streamline Length	0.5 m

Qui Resolution controlla il numero di intervalli dei semi lungo la linea dei semi; non è la risoluzione della mesh CFD.

In questo caso stazionario le linee di corrente coincidono con le traiettorie delle particelle. Devono essere rette e parallele alla direzione x.

19. Risoluzione dei problemi generici

Sintomo	Controllo / correzione
blockMesh o checkMesh non riescono ad avviare il caso	verificare che esista system/controlDict
foamPostProcess non crea la cartella postProcessing	richiamare esplicitamente la funzione con -func graphUniform
il file atteso line_U.xy non è presente	OpenFOAM 14 scrive line.xy con questa configurazione
l’include di graphUniform non viene trovato	usare caseDicts/postProcessing/graphs/graphUniform.cfg
compare solo un piccolo numero di linee di corrente	verificare che l’intera linea dei semi stia dentro il dominio CFD
la vista delle linee di corrente è visivamente confusa	visualizzare la sorgente CFD come Outline o ridurre l’opacità della superficie

20. Criteri di accettazione

Criterio di accettazione	Stato
Ambiente OpenFOAM operativo	SUPERATO
Mesh generata e checkMesh riuscito	SUPERATO
Il calcolo stazionario arriva a un End normale	SUPERATO
Profilo di velocità simmetrico	SUPERATO
Uy trascurabile e Uz = 0	SUPERATO
RMSE normalizzato < 0,2%	SUPERATO – 0,150343%
La pressione diminuisce lungo il canale	SUPERATO
Visualizzazione ParaView client-server operativa	SUPERATO
Le linee di corrente attraversano il canale e restano parallele a x	SUPERATO

Stato finale: CASO 01 COMPLETATO E VALIDATO.

21. Checklist di riproducibilità

• Caricare l’ambiente di OpenFOAM 14.
• Creare $FOAM_RUN/01_laminar_channel.
• Creare blockMeshDict e controlDict.
• Eseguire blockMesh e checkMesh.
• Creare 0/U e 0/p.
• Creare constant/physicalProperties e constant/momentumTransport.
• Creare fvSchemes e fvSolution.
• Eseguire foamRun | tee log.foamRun.
• Campionare U a x = 0,45 m con graphUniform.
• Confrontare il profilo campionato con la teoria di Poiseuille.
• Facoltativamente, campionare p lungo x e confrontarla con il gradiente di pressione analitico.
• Creare il file marcatore .foam.
• Avviare pvserver e connettersi da ParaView 5.13.2.
• Esaminare U, p e le linee di corrente all’etichetta di iterazione finale.

Questo caso è il riferimento validato per le simulazioni CFD successive del percorso: nel caso 02 lo stesso canale serve per lo studio di convergenza della mesh.