diff --git a/.gitignore b/.gitignore index 77ce0d31..e1e609d4 100644 --- a/.gitignore +++ b/.gitignore @@ -128,5 +128,9 @@ dmypy.json # Pyre type checker .pyre/ +.viminfo +.bash_history + # output .DS_Store +*.stl diff --git a/OFsolvers/tutorial_cases/loop_reactor/0.orig/T.air b/OFsolvers/tutorial_cases/loop_reactor/0.orig/T.air new file mode 100644 index 00000000..b71958b1 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/0.orig/T.air @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class volScalarField; + object T.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 1 0 0 0]; + +internalField uniform 300; + +boundaryField +{ + outlet + { + type inletOutlet; + phi phi.air; + inletValue $internalField; + value $internalField; + } + inlet + { + type fixedValue; + value $internalField; + } + defaultFaces + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/0.orig/T.water b/OFsolvers/tutorial_cases/loop_reactor/0.orig/T.water new file mode 100644 index 00000000..da0d6ff6 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/0.orig/T.water @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class volScalarField; + object T.water; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 1 0 0 0]; + +internalField uniform 350; + +boundaryField +{ + outlet + { + type inletOutlet; + phi phi.water; + inletValue $internalField; + value $internalField; + } + inlet + { + type fixedValue; + value $internalField; + } + defaultFaces + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/0.orig/U.air b/OFsolvers/tutorial_cases/loop_reactor/0.orig/U.air new file mode 100644 index 00000000..8bf67bd0 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/0.orig/U.air @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format binary; + class volVectorField; + object U.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0.0 0); + +boundaryField +{ + inlet + { + type fixedValue; + value uniform (0.0 0.025 0.0); + } + outlet + { + type pressureInletOutletVelocity; + phi phi.air; + value $internalField; + } + defaultFaces + { + type fixedValue; + value uniform (0 0 0); + } +} + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/0.orig/U.water b/OFsolvers/tutorial_cases/loop_reactor/0.orig/U.water new file mode 100644 index 00000000..de2877d1 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/0.orig/U.water @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format binary; + class volVectorField; + object U.water; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + inlet + { + type fixedValue; + value $internalField; + } + outlet + { + type pressureInletOutletVelocity; + phi phi.water; + value $internalField; + } + defaultFaces + { + type fixedValue; + value uniform (0 0 0); + } +} + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/0.orig/alpha.air.orig b/OFsolvers/tutorial_cases/loop_reactor/0.orig/alpha.air.orig new file mode 100644 index 00000000..576f080f --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/0.orig/alpha.air.orig @@ -0,0 +1,41 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class volScalarField; + location "0"; + object alpha.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + inlet + { + type fixedValue; + value uniform 0.5; + } + outlet + { + type inletOutlet; + phi phi.air; + inletValue uniform 1; + value uniform 1; + } + defaultFaces + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/0.orig/alpha.water.orig b/OFsolvers/tutorial_cases/loop_reactor/0.orig/alpha.water.orig new file mode 100644 index 00000000..b39ce2b9 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/0.orig/alpha.water.orig @@ -0,0 +1,40 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class volScalarField; + object alpha.water; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 0 0 0 0 0 0]; + +internalField uniform 1; + +boundaryField +{ + inlet + { + type fixedValue; + value uniform 0.5; + } + outlet + { + type inletOutlet; + phi phi.water; + inletValue uniform 0; + value uniform 0; + } + defaultFaces + { + type zeroGradient; + } +} + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/0.orig/p b/OFsolvers/tutorial_cases/loop_reactor/0.orig/p new file mode 100644 index 00000000..b351061b --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/0.orig/p @@ -0,0 +1,39 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class volScalarField; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -1 -2 0 0 0 0]; + +internalField uniform 1e5; + +boundaryField +{ + inlet + { + type calculated; + value $internalField; + } + outlet + { + type calculated; + value $internalField; + } + defaultFaces + { + type calculated; + value $internalField; + } +} + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/0.orig/p_rgh b/OFsolvers/tutorial_cases/loop_reactor/0.orig/p_rgh new file mode 100644 index 00000000..04d06df4 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/0.orig/p_rgh @@ -0,0 +1,43 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class volScalarField; + object p_rgh; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [1 -1 -2 0 0 0 0]; + +internalField uniform 1e5; + +boundaryField +{ + inlet + { + type fixedFluxPressure; + value $internalField; + } + outlet + { + type prghTotalPressure; + p0 $internalField; + U U.air; + phi phi.air; + rho thermo:rho.air; + value $internalField; + } + defaultFaces + { + type fixedFluxPressure; + value $internalField; + } +} + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/Allclean b/OFsolvers/tutorial_cases/loop_reactor/Allclean new file mode 100755 index 00000000..4003a12b --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/Allclean @@ -0,0 +1,15 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory + +# Source tutorial clean functions +. $WM_PROJECT_DIR/bin/tools/CleanFunctions + +# Remove surface, features and solution +#rm -rf constant/extendedFeatureEdgeMesh > /dev/null 2>&1 +#rm -f constant/triSurface/*.eMesh > /dev/null 2>&1 +#rm -rf constant/polyMesh > /dev/null 2>&1 +#rm -rf processor* > /dev/null 2>&1 +rm -rf 0 +cleanCase + +#------------------------------------------------------------------------------ diff --git a/OFsolvers/tutorial_cases/loop_reactor/constant/g b/OFsolvers/tutorial_cases/loop_reactor/constant/g new file mode 100644 index 00000000..770a5619 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/constant/g @@ -0,0 +1,21 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class uniformDimensionedVectorField; + location "constant"; + object g; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -2 0 0 0 0]; +value (0 -9.81 0); + + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/constant/momentumTransport.air b/OFsolvers/tutorial_cases/loop_reactor/constant/momentumTransport.air new file mode 100644 index 00000000..bceba9af --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/constant/momentumTransport.air @@ -0,0 +1,19 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "constant"; + object momentumTransport.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType laminar; + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/constant/momentumTransport.water b/OFsolvers/tutorial_cases/loop_reactor/constant/momentumTransport.water new file mode 100644 index 00000000..05a1ed6f --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/constant/momentumTransport.water @@ -0,0 +1,19 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "constant"; + object momentumTransport.water; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType laminar; + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/constant/phaseProperties b/OFsolvers/tutorial_cases/loop_reactor/constant/phaseProperties new file mode 100644 index 00000000..ba8e290d --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/constant/phaseProperties @@ -0,0 +1,147 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "constant"; + object phaseProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +type basicMultiphaseSystem; + +phases (air water); + +referencePhase water; + +air +{ + type purePhaseModel; + diameterModel isothermal; + isothermalCoeffs + { + d0 3e-3; + p0 1e5; + } + + residualAlpha 1e-6; +} + +water +{ + type purePhaseModel; + diameterModel constant; + constantCoeffs + { + d 1e-4; + } + + residualAlpha 1e-6; +} + +blending +{ + default + { + type linear; + minFullyContinuousAlpha.air 0.7; + minPartlyContinuousAlpha.air 0.3; + minFullyContinuousAlpha.water 0.7; + minPartlyContinuousAlpha.water 0.3; + } + + drag + { + type linear; + minFullyContinuousAlpha.air 0.7; + minPartlyContinuousAlpha.air 0.5; + minFullyContinuousAlpha.water 0.7; + minPartlyContinuousAlpha.water 0.5; + } +} + +surfaceTension +( + (air and water) + { + type constant; + sigma 0.07; + } +); + +interfaceCompression +(); + +aspectRatio +( + (air in water) + { + type Wellek; + } +); + +drag +( + (air in water) + { + type SchillerNaumann; + residualRe 1e-3; + } + + (water in air) + { + type SchillerNaumann; + residualRe 1e-3; + } + +); + +virtualMass +( + (air in water) + { + type constantCoefficient; + Cvm 0.5; + } + + (water in air) + { + type constantCoefficient; + Cvm 0.5; + } +); + +heatTransfer +( + (air in water) + { + type RanzMarshall; + residualAlpha 1e-4; + } + + (water in air) + { + type RanzMarshall; + residualAlpha 1e-4; + } +); + +phaseTransfer +(); + +lift +(); + +wallLubrication +(); + +turbulentDispersion +(); + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/constant/thermophysicalProperties.air b/OFsolvers/tutorial_cases/loop_reactor/constant/thermophysicalProperties.air new file mode 100644 index 00000000..8b6595ef --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/constant/thermophysicalProperties.air @@ -0,0 +1,47 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "constant"; + object physicalProperties.air; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType +{ + type heRhoThermo; + mixture pureMixture; + transport const; + thermo hConst; + equationOfState perfectGas; + specie specie; + energy sensibleInternalEnergy; +} + +mixture +{ + specie + { + molWeight 28.9; + } + thermodynamics + { + Cp 1007; + Hf 0; + } + transport + { + mu 1.84e-05; + Pr 0.7; + } +} + + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/constant/thermophysicalProperties.water b/OFsolvers/tutorial_cases/loop_reactor/constant/thermophysicalProperties.water new file mode 100644 index 00000000..8843e086 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/constant/thermophysicalProperties.water @@ -0,0 +1,51 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "constant"; + object physicalProperties.water; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +thermoType +{ + type heRhoThermo; + mixture pureMixture; + transport const; + thermo eConst; + equationOfState rPolynomial; + specie specie; + energy sensibleInternalEnergy; +} + +mixture +{ + specie + { + molWeight 18; + } + equationOfState + { + C (0.001278 -2.1055e-06 3.9689e-09 4.3772e-13 -2.0225e-16); + } + thermodynamics + { + Cv 4195; + Hf 0; + } + transport + { + mu 3.645e-4; + Pr 2.289; + } +} + + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/run.sh b/OFsolvers/tutorial_cases/loop_reactor/run.sh new file mode 100644 index 00000000..e9f4b557 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/run.sh @@ -0,0 +1,38 @@ +# Generate blockmeshDict +python ../../../applications/write_block_rect_mesh.py -i ../../../bird/meshing/block_rect_mesh_templates/loopReactor/input.json -o system + +# Generate boundary stl +python ../../../applications/write_stl_patch.py -i ../../../bird/preProcess/stl_patch/bc_patch_mesh_template/loop_reactor/inlets_outlets.json + +# Clean case +./Allclean + +# Mesh gen +blockMesh -dict system/blockMeshDict + +# Inlet BC +surfaceToPatch -tol 1e-3 inlets.stl +export newmeshdir=$(foamListTimes -latestTime) +rm -rf constant/polyMesh/ +cp -r $newmeshdir/polyMesh ./constant +rm -rf $newmeshdir +sed -i 's/inlets\.stl/inlet/g' ./constant/polyMesh/boundary + +# Outlet BC +surfaceToPatch -tol 1e-3 outlets.stl +export newmeshdir=$(foamListTimes -latestTime) +rm -rf constant/polyMesh/ +cp -r $newmeshdir/polyMesh ./constant +rm -rf $newmeshdir +sed -i 's/outlets\.stl/outlet/g' ./constant/polyMesh/boundary + + +# Scale +transformPoints "scale=(0.05 0.05 0.05)" + +# setup IC +cp -r 0.orig 0 +setFields + +# Run +multiphaseEulerFoam diff --git a/OFsolvers/tutorial_cases/loop_reactor/system/blockMeshDict b/OFsolvers/tutorial_cases/loop_reactor/system/blockMeshDict new file mode 100644 index 00000000..94a55474 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/system/blockMeshDict @@ -0,0 +1,594 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} + +convertToMeters 1.0; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +vertices +( +( 0.0 0.0 0.0) +( 1.0 0.0 0.0) +( 2.0 0.0 0.0) +( 3.0 0.0 0.0) +( 4.0 0.0 0.0) +( 5.0 0.0 0.0) +( 6.0 0.0 0.0) +( 7.0 0.0 0.0) +( 8.0 0.0 0.0) +( 9.0 0.0 0.0) +( 10.0 0.0 0.0) +( 0.0 1.0 0.0) +( 1.0 1.0 0.0) +( 2.0 1.0 0.0) +( 3.0 1.0 0.0) +( 4.0 1.0 0.0) +( 5.0 1.0 0.0) +( 6.0 1.0 0.0) +( 7.0 1.0 0.0) +( 8.0 1.0 0.0) +( 9.0 1.0 0.0) +( 10.0 1.0 0.0) +( 0.0 2.0 0.0) +( 1.0 2.0 0.0) +( 2.0 2.0 0.0) +( 3.0 2.0 0.0) +( 4.0 2.0 0.0) +( 5.0 2.0 0.0) +( 6.0 2.0 0.0) +( 7.0 2.0 0.0) +( 8.0 2.0 0.0) +( 9.0 2.0 0.0) +( 10.0 2.0 0.0) +( 0.0 3.0 0.0) +( 1.0 3.0 0.0) +( 2.0 3.0 0.0) +( 3.0 3.0 0.0) +( 4.0 3.0 0.0) +( 5.0 3.0 0.0) +( 6.0 3.0 0.0) +( 7.0 3.0 0.0) +( 8.0 3.0 0.0) +( 9.0 3.0 0.0) +( 10.0 3.0 0.0) +( 0.0 4.0 0.0) +( 1.0 4.0 0.0) +( 2.0 4.0 0.0) +( 3.0 4.0 0.0) +( 4.0 4.0 0.0) +( 5.0 4.0 0.0) +( 6.0 4.0 0.0) +( 7.0 4.0 0.0) +( 8.0 4.0 0.0) +( 9.0 4.0 0.0) +( 10.0 4.0 0.0) +( 0.0 5.0 0.0) +( 1.0 5.0 0.0) +( 2.0 5.0 0.0) +( 3.0 5.0 0.0) +( 4.0 5.0 0.0) +( 5.0 5.0 0.0) +( 6.0 5.0 0.0) +( 7.0 5.0 0.0) +( 8.0 5.0 0.0) +( 9.0 5.0 0.0) +( 10.0 5.0 0.0) +( 0.0 0.0 1.0) +( 1.0 0.0 1.0) +( 2.0 0.0 1.0) +( 3.0 0.0 1.0) +( 4.0 0.0 1.0) +( 5.0 0.0 1.0) +( 6.0 0.0 1.0) +( 7.0 0.0 1.0) +( 8.0 0.0 1.0) +( 9.0 0.0 1.0) +( 10.0 0.0 1.0) +( 0.0 1.0 1.0) +( 1.0 1.0 1.0) +( 2.0 1.0 1.0) +( 3.0 1.0 1.0) +( 4.0 1.0 1.0) +( 5.0 1.0 1.0) +( 6.0 1.0 1.0) +( 7.0 1.0 1.0) +( 8.0 1.0 1.0) +( 9.0 1.0 1.0) +( 10.0 1.0 1.0) +( 0.0 2.0 1.0) +( 1.0 2.0 1.0) +( 2.0 2.0 1.0) +( 3.0 2.0 1.0) +( 4.0 2.0 1.0) +( 5.0 2.0 1.0) +( 6.0 2.0 1.0) +( 7.0 2.0 1.0) +( 8.0 2.0 1.0) +( 9.0 2.0 1.0) +( 10.0 2.0 1.0) +( 0.0 3.0 1.0) +( 1.0 3.0 1.0) +( 2.0 3.0 1.0) +( 3.0 3.0 1.0) +( 4.0 3.0 1.0) +( 5.0 3.0 1.0) +( 6.0 3.0 1.0) +( 7.0 3.0 1.0) +( 8.0 3.0 1.0) +( 9.0 3.0 1.0) +( 10.0 3.0 1.0) +( 0.0 4.0 1.0) +( 1.0 4.0 1.0) +( 2.0 4.0 1.0) +( 3.0 4.0 1.0) +( 4.0 4.0 1.0) +( 5.0 4.0 1.0) +( 6.0 4.0 1.0) +( 7.0 4.0 1.0) +( 8.0 4.0 1.0) +( 9.0 4.0 1.0) +( 10.0 4.0 1.0) +( 0.0 5.0 1.0) +( 1.0 5.0 1.0) +( 2.0 5.0 1.0) +( 3.0 5.0 1.0) +( 4.0 5.0 1.0) +( 5.0 5.0 1.0) +( 6.0 5.0 1.0) +( 7.0 5.0 1.0) +( 8.0 5.0 1.0) +( 9.0 5.0 1.0) +( 10.0 5.0 1.0) +( 0.0 0.0 2.0) +( 1.0 0.0 2.0) +( 2.0 0.0 2.0) +( 3.0 0.0 2.0) +( 4.0 0.0 2.0) +( 5.0 0.0 2.0) +( 6.0 0.0 2.0) +( 7.0 0.0 2.0) +( 8.0 0.0 2.0) +( 9.0 0.0 2.0) +( 10.0 0.0 2.0) +( 0.0 1.0 2.0) +( 1.0 1.0 2.0) +( 2.0 1.0 2.0) +( 3.0 1.0 2.0) +( 4.0 1.0 2.0) +( 5.0 1.0 2.0) +( 6.0 1.0 2.0) +( 7.0 1.0 2.0) +( 8.0 1.0 2.0) +( 9.0 1.0 2.0) +( 10.0 1.0 2.0) +( 0.0 2.0 2.0) +( 1.0 2.0 2.0) +( 2.0 2.0 2.0) +( 3.0 2.0 2.0) +( 4.0 2.0 2.0) +( 5.0 2.0 2.0) +( 6.0 2.0 2.0) +( 7.0 2.0 2.0) +( 8.0 2.0 2.0) +( 9.0 2.0 2.0) +( 10.0 2.0 2.0) +( 0.0 3.0 2.0) +( 1.0 3.0 2.0) +( 2.0 3.0 2.0) +( 3.0 3.0 2.0) +( 4.0 3.0 2.0) +( 5.0 3.0 2.0) +( 6.0 3.0 2.0) +( 7.0 3.0 2.0) +( 8.0 3.0 2.0) +( 9.0 3.0 2.0) +( 10.0 3.0 2.0) +( 0.0 4.0 2.0) +( 1.0 4.0 2.0) +( 2.0 4.0 2.0) +( 3.0 4.0 2.0) +( 4.0 4.0 2.0) +( 5.0 4.0 2.0) +( 6.0 4.0 2.0) +( 7.0 4.0 2.0) +( 8.0 4.0 2.0) +( 9.0 4.0 2.0) +( 10.0 4.0 2.0) +( 0.0 5.0 2.0) +( 1.0 5.0 2.0) +( 2.0 5.0 2.0) +( 3.0 5.0 2.0) +( 4.0 5.0 2.0) +( 5.0 5.0 2.0) +( 6.0 5.0 2.0) +( 7.0 5.0 2.0) +( 8.0 5.0 2.0) +( 9.0 5.0 2.0) +( 10.0 5.0 2.0) +( 0.0 0.0 3.0) +( 1.0 0.0 3.0) +( 2.0 0.0 3.0) +( 3.0 0.0 3.0) +( 4.0 0.0 3.0) +( 5.0 0.0 3.0) +( 6.0 0.0 3.0) +( 7.0 0.0 3.0) +( 8.0 0.0 3.0) +( 9.0 0.0 3.0) +( 10.0 0.0 3.0) +( 0.0 1.0 3.0) +( 1.0 1.0 3.0) +( 2.0 1.0 3.0) +( 3.0 1.0 3.0) +( 4.0 1.0 3.0) +( 5.0 1.0 3.0) +( 6.0 1.0 3.0) +( 7.0 1.0 3.0) +( 8.0 1.0 3.0) +( 9.0 1.0 3.0) +( 10.0 1.0 3.0) +( 0.0 2.0 3.0) +( 1.0 2.0 3.0) +( 2.0 2.0 3.0) +( 3.0 2.0 3.0) +( 4.0 2.0 3.0) +( 5.0 2.0 3.0) +( 6.0 2.0 3.0) +( 7.0 2.0 3.0) +( 8.0 2.0 3.0) +( 9.0 2.0 3.0) +( 10.0 2.0 3.0) +( 0.0 3.0 3.0) +( 1.0 3.0 3.0) +( 2.0 3.0 3.0) +( 3.0 3.0 3.0) +( 4.0 3.0 3.0) +( 5.0 3.0 3.0) +( 6.0 3.0 3.0) +( 7.0 3.0 3.0) +( 8.0 3.0 3.0) +( 9.0 3.0 3.0) +( 10.0 3.0 3.0) +( 0.0 4.0 3.0) +( 1.0 4.0 3.0) +( 2.0 4.0 3.0) +( 3.0 4.0 3.0) +( 4.0 4.0 3.0) +( 5.0 4.0 3.0) +( 6.0 4.0 3.0) +( 7.0 4.0 3.0) +( 8.0 4.0 3.0) +( 9.0 4.0 3.0) +( 10.0 4.0 3.0) +( 0.0 5.0 3.0) +( 1.0 5.0 3.0) +( 2.0 5.0 3.0) +( 3.0 5.0 3.0) +( 4.0 5.0 3.0) +( 5.0 5.0 3.0) +( 6.0 5.0 3.0) +( 7.0 5.0 3.0) +( 8.0 5.0 3.0) +( 9.0 5.0 3.0) +( 10.0 5.0 3.0) +( 0.0 0.0 4.0) +( 1.0 0.0 4.0) +( 2.0 0.0 4.0) +( 3.0 0.0 4.0) +( 4.0 0.0 4.0) +( 5.0 0.0 4.0) +( 6.0 0.0 4.0) +( 7.0 0.0 4.0) +( 8.0 0.0 4.0) +( 9.0 0.0 4.0) +( 10.0 0.0 4.0) +( 0.0 1.0 4.0) +( 1.0 1.0 4.0) +( 2.0 1.0 4.0) +( 3.0 1.0 4.0) +( 4.0 1.0 4.0) +( 5.0 1.0 4.0) +( 6.0 1.0 4.0) +( 7.0 1.0 4.0) +( 8.0 1.0 4.0) +( 9.0 1.0 4.0) +( 10.0 1.0 4.0) +( 0.0 2.0 4.0) +( 1.0 2.0 4.0) +( 2.0 2.0 4.0) +( 3.0 2.0 4.0) +( 4.0 2.0 4.0) +( 5.0 2.0 4.0) +( 6.0 2.0 4.0) +( 7.0 2.0 4.0) +( 8.0 2.0 4.0) +( 9.0 2.0 4.0) +( 10.0 2.0 4.0) +( 0.0 3.0 4.0) +( 1.0 3.0 4.0) +( 2.0 3.0 4.0) +( 3.0 3.0 4.0) +( 4.0 3.0 4.0) +( 5.0 3.0 4.0) +( 6.0 3.0 4.0) +( 7.0 3.0 4.0) +( 8.0 3.0 4.0) +( 9.0 3.0 4.0) +( 10.0 3.0 4.0) +( 0.0 4.0 4.0) +( 1.0 4.0 4.0) +( 2.0 4.0 4.0) +( 3.0 4.0 4.0) +( 4.0 4.0 4.0) +( 5.0 4.0 4.0) +( 6.0 4.0 4.0) +( 7.0 4.0 4.0) +( 8.0 4.0 4.0) +( 9.0 4.0 4.0) +( 10.0 4.0 4.0) +( 0.0 5.0 4.0) +( 1.0 5.0 4.0) +( 2.0 5.0 4.0) +( 3.0 5.0 4.0) +( 4.0 5.0 4.0) +( 5.0 5.0 4.0) +( 6.0 5.0 4.0) +( 7.0 5.0 4.0) +( 8.0 5.0 4.0) +( 9.0 5.0 4.0) +( 10.0 5.0 4.0) +( 0.0 0.0 5.0) +( 1.0 0.0 5.0) +( 2.0 0.0 5.0) +( 3.0 0.0 5.0) +( 4.0 0.0 5.0) +( 5.0 0.0 5.0) +( 6.0 0.0 5.0) +( 7.0 0.0 5.0) +( 8.0 0.0 5.0) +( 9.0 0.0 5.0) +( 10.0 0.0 5.0) +( 0.0 1.0 5.0) +( 1.0 1.0 5.0) +( 2.0 1.0 5.0) +( 3.0 1.0 5.0) +( 4.0 1.0 5.0) +( 5.0 1.0 5.0) +( 6.0 1.0 5.0) +( 7.0 1.0 5.0) +( 8.0 1.0 5.0) +( 9.0 1.0 5.0) +( 10.0 1.0 5.0) +( 0.0 2.0 5.0) +( 1.0 2.0 5.0) +( 2.0 2.0 5.0) +( 3.0 2.0 5.0) +( 4.0 2.0 5.0) +( 5.0 2.0 5.0) +( 6.0 2.0 5.0) +( 7.0 2.0 5.0) +( 8.0 2.0 5.0) +( 9.0 2.0 5.0) +( 10.0 2.0 5.0) +( 0.0 3.0 5.0) +( 1.0 3.0 5.0) +( 2.0 3.0 5.0) +( 3.0 3.0 5.0) +( 4.0 3.0 5.0) +( 5.0 3.0 5.0) +( 6.0 3.0 5.0) +( 7.0 3.0 5.0) +( 8.0 3.0 5.0) +( 9.0 3.0 5.0) +( 10.0 3.0 5.0) +( 0.0 4.0 5.0) +( 1.0 4.0 5.0) +( 2.0 4.0 5.0) +( 3.0 4.0 5.0) +( 4.0 4.0 5.0) +( 5.0 4.0 5.0) +( 6.0 4.0 5.0) +( 7.0 4.0 5.0) +( 8.0 4.0 5.0) +( 9.0 4.0 5.0) +( 10.0 4.0 5.0) +( 0.0 5.0 5.0) +( 1.0 5.0 5.0) +( 2.0 5.0 5.0) +( 3.0 5.0 5.0) +( 4.0 5.0 5.0) +( 5.0 5.0 5.0) +( 6.0 5.0 5.0) +( 7.0 5.0 5.0) +( 8.0 5.0 5.0) +( 9.0 5.0 5.0) +( 10.0 5.0 5.0) +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +blocks +( + + //block 0 +hex (0 1 12 11 66 67 78 77 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 1 +hex (1 2 13 12 67 68 79 78 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 2 +hex (2 3 14 13 68 69 80 79 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 3 +hex (3 4 15 14 69 70 81 80 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 4 +hex (4 5 16 15 70 71 82 81 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 5 +hex (5 6 17 16 71 72 83 82 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 6 +hex (6 7 18 17 72 73 84 83 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 7 +hex (7 8 19 18 73 74 85 84 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 8 +hex (8 9 20 19 74 75 86 85 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 9 +hex (9 10 21 20 75 76 87 86 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 10 +hex (75 76 87 86 141 142 153 152 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 11 +hex (141 142 153 152 207 208 219 218 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 12 +hex (207 208 219 218 273 274 285 284 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 13 +hex (273 274 285 284 339 340 351 350 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 14 +hex (272 273 284 283 338 339 350 349 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 15 +hex (271 272 283 282 337 338 349 348 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 16 +hex (270 271 282 281 336 337 348 347 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 17 +hex (269 270 281 280 335 336 347 346 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 18 +hex (268 269 280 279 334 335 346 345 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 19 +hex (267 268 279 278 333 334 345 344 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 20 +hex (266 267 278 277 332 333 344 343 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 21 +hex (265 266 277 276 331 332 343 342 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 22 +hex (264 265 276 275 330 331 342 341 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 23 +hex (275 276 287 286 341 342 353 352 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 24 +hex (286 287 298 297 352 353 364 363 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 25 +hex (297 298 309 308 363 364 375 374 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 26 +hex (308 309 320 319 374 375 386 385 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 27 +hex (242 243 254 253 308 309 320 319 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 28 +hex (176 177 188 187 242 243 254 253 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 29 +hex (110 111 122 121 176 177 188 187 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 30 +hex (44 45 56 55 110 111 122 121 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 31 +hex (33 34 45 44 99 100 111 110 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 32 +hex (22 23 34 33 88 89 100 99 ) +( 10 10 10 ) +SimpleGrading (1 1 1) + + //block 33 +hex (11 12 23 22 77 78 89 88 ) +( 10 10 10 ) +SimpleGrading (1 1 1) +); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +defaultPatch +{ type wall;} + +patches +( +); diff --git a/OFsolvers/tutorial_cases/loop_reactor/system/controlDict b/OFsolvers/tutorial_cases/loop_reactor/system/controlDict new file mode 100644 index 00000000..66306ff7 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/system/controlDict @@ -0,0 +1,59 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application multiphaseEulerFoam; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 100; + +deltaT 0.005; + +writeControl runTime; + +writeInterval 1; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression off; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable yes; + +adjustTimeStep yes; + +maxCo 0.5; + +maxDeltaT 1; + +//functions +//{ +// #includeFunc fieldAverage(U.air, U.water, alpha.air, p) +//} + + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/system/decomposeParDict b/OFsolvers/tutorial_cases/loop_reactor/system/decomposeParDict new file mode 100755 index 00000000..a1cfef65 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/system/decomposeParDict @@ -0,0 +1,30 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: 3.0.x | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object decomposeParDict; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 16; + +method hierarchical; + +hierarchicalCoeffs +{ + n (4 4 1); + delta 0.001; + order xyz; +} + + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/system/fvConstraints b/OFsolvers/tutorial_cases/loop_reactor/system/fvConstraints new file mode 100644 index 00000000..164c3e0b --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/system/fvConstraints @@ -0,0 +1,23 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + object fvConstraints; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +limitp +{ + type limitPressure; + + min 1e4; +} + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/system/fvSchemes b/OFsolvers/tutorial_cases/loop_reactor/system/fvSchemes new file mode 100644 index 00000000..806a9d3d --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/system/fvSchemes @@ -0,0 +1,63 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{ + default Euler; +} + +gradSchemes +{ + default Gauss linear; +} + +divSchemes +{ + default none; + + "div\(phi,alpha.*\)" Gauss vanLeer; + + "div\(phir,alpha.*,alpha.*\)" Gauss vanLeer; + + "div\(alphaRhoPhi.*,U.*\)" Gauss limitedLinearV 1; + "div\(phi.*,U.*\)" Gauss limitedLinearV 1; + + "div\(alphaRhoPhi.*,e.*\)" Gauss limitedLinear 1; + "div\(alphaRhoPhi.*,K.*\)" Gauss limitedLinear 1; + "div\(alphaRhoPhi.*,\(p\|thermo:rho.*\)\)" Gauss limitedLinear 1; + + + + "div\(\(\(\(alpha.*\*thermo:rho.*\)*nuEff.*\)*dev2\(T\(grad\(U.*\)\)\)\)\)" Gauss linear; +} + +laplacianSchemes +{ + default Gauss linear uncorrected; +} + +interpolationSchemes +{ + default linear; +} + +snGradSchemes +{ + default uncorrected; +} + + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/system/fvSolution b/OFsolvers/tutorial_cases/loop_reactor/system/fvSolution new file mode 100644 index 00000000..9523801f --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/system/fvSolution @@ -0,0 +1,75 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + "alpha.*" + { + nAlphaCorr 1; + nAlphaSubCycles 2; + } + + p_rgh + { + solver GAMG; + smoother DIC; + tolerance 1e-8; + relTol 0; + } + + p_rghFinal + { + $p_rgh; + relTol 0; + } + + "U.*" + { + solver smoothSolver; + smoother symGaussSeidel; + tolerance 1e-5; + relTol 0; + minIter 1; + } + + "e.*" + { + solver smoothSolver; + smoother symGaussSeidel; + tolerance 1e-8; + relTol 0; + minIter 1; + } +} + +PIMPLE +{ + nOuterCorrectors 3; + nCorrectors 1; + nNonOrthogonalCorrectors 0; + +} + +relaxationFactors +{ + equations + { + ".*" 1; + } +} + + +// ************************************************************************* // diff --git a/OFsolvers/tutorial_cases/loop_reactor/system/setFieldsDict b/OFsolvers/tutorial_cases/loop_reactor/system/setFieldsDict new file mode 100644 index 00000000..97350461 --- /dev/null +++ b/OFsolvers/tutorial_cases/loop_reactor/system/setFieldsDict @@ -0,0 +1,37 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: dev + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + format ascii; + class dictionary; + location "system"; + object setFieldsDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +defaultFieldValues +( + volScalarFieldValue alpha.air 1 + volScalarFieldValue alpha.water 0 +); + +regions +( + boxToCell + { + box (-1.0 -1.0 -1.0) (10 0.15 10.0); + fieldValues + ( + volScalarFieldValue alpha.air 0 + volScalarFieldValue alpha.water 1 + ); + } +); + + +// ************************************************************************* // diff --git a/README.md b/README.md index 2473942f..f15aba1f 100644 --- a/README.md +++ b/README.md @@ -53,35 +53,6 @@ Mesh visualized in Paraview

-### Generate STL mesh - -`python applications/write_stl_mesh.py -v -cr 0.25 -na 12 -aw 0.1 -al 0.5` - -Generates - -

- -

- - -### Manual - -``` -usage: write_stl_mesh.py [-h] [-cr] [-na] [-aw] [-al] [-v] - -Generate Spider Sparger STL - -optional arguments: - -h, --help show this help message and exit - -cr , --centerRadius - Radius of the center distributor - -na , --nArms Number of spider arms - -aw , --armsWidth Width of spider arms - -al , --armsLength Length of spider arms - -v, --verbose plot on screen - -``` - ### Block cylindrical meshing Generates `blockMeshDict` in `system` @@ -240,13 +211,57 @@ The corners are defined in the `input.json` } ... ``` +To see how to use this on an actual case see `OFsolvers/tutorial_cases/loop_reactor` + +## Preprocess +### Generate STL mesh + +Boundaries may be specified with `surfaceToPatch` utility in OpenFOAM, based on STL files that can be generated with + +`python applications/write_stl_patch.py -v` + +Generates + +

+ +

+ + +To see how to use this on an actual case see `OFsolvers/tutorial_cases/loop_reactor` + +### Manual + +``` +usage: write_stl_patch.py [-h] [-i] [-v] + +Generate boundary patch + +options: + -h, --help show this help message and exit + -i , --input Boundary patch Json input + -v, --verbose plot on screen +``` + +### How to change the set of shapes in the boundary patch +Edit the json files read when generating the mesh. In the case below, the boundary condition `inlets` consists of 3 discs + +``` +{ + "inlets": [ + {"type": "circle", "centx": 5.0, "centy": 0.0, "centz": 0.5, "radius": 0.4, "normal_dir": 1,"nelements": 50}, + {"type": "circle", "centx": 2.5, "centy": 0.0, "centz": 0.5, "radius": 0.4, "normal_dir": 1,"nelements": 50}, + {"type": "circle", "centx": 7.5, "centy": 0.0, "centz": 0.5, "radius": 0.4, "normal_dir": 1,"nelements": 50} + ], +} +... +``` ## Postprocess ### Perform early prediction -`python applications/earlyPredicition.py -df bird/postProcess/data_early` +`python applications/early_prediction.py -df bird/postProcess/data_early` Generates diff --git a/applications/compute_conditional_mean.py b/applications/compute_conditional_mean.py index 27a1517d..cab87e25 100644 --- a/applications/compute_conditional_mean.py +++ b/applications/compute_conditional_mean.py @@ -10,7 +10,6 @@ def main(): - parser = argparse.ArgumentParser( description="Compute conditional means of OpenFOAM fields" ) diff --git a/applications/write_stl_mesh.py b/applications/write_stl_mesh.py deleted file mode 100644 index 45a18b26..00000000 --- a/applications/write_stl_mesh.py +++ /dev/null @@ -1,68 +0,0 @@ -import argparse -import sys - -import numpy as np - -from bird.meshing.stl_mesh_tools import makeSpider, saveSTL - -if __name__ == "__main__": - parser = argparse.ArgumentParser(description="Generate Spider Sparger STL") - parser.add_argument( - "-cr", - "--centerRadius", - type=float, - metavar="", - required=False, - help="Radius of the center distributor", - default=0.25, - ) - parser.add_argument( - "-na", - "--nArms", - type=int, - metavar="", - required=False, - help="Number of spider arms", - default=12, - ) - parser.add_argument( - "-aw", - "--armsWidth", - type=float, - metavar="", - required=False, - help="Width of spider arms", - default=0.1, - ) - parser.add_argument( - "-al", - "--armsLength", - type=float, - metavar="", - required=False, - help="Length of spider arms", - default=0.5, - ) - parser.add_argument( - "-v", "--verbose", action="store_true", help="plot on screen" - ) - - args = parser.parse_args() - # Spider - combined, globalArea = makeSpider( - centerRad=args.centerRadius, - nArms=args.nArms, - widthArms=args.armsWidth, - lengthArms=args.armsLength, - ) - print(f"\tglobalArea = {globalArea}") - - saveSTL(combined, filename="spg.stl") - - if args.verbose: - # plot - from bird.utilities.stl_plotting import plotSTL, plt, pretty_labels - - axes = plotSTL("spg.stl") - pretty_labels("x", "y", zlabel="z", fontsize=14, ax=axes) - plt.show() diff --git a/applications/write_stl_patch.py b/applications/write_stl_patch.py new file mode 100644 index 00000000..3956e4ed --- /dev/null +++ b/applications/write_stl_patch.py @@ -0,0 +1,39 @@ +import os + +import numpy as np +import stl + +from bird import BIRD_PRE_PATCH_TEMP_DIR +from bird.meshing._mesh_tools import parseJsonFile +from bird.preProcess.stl_patch.stl_bc import write_boundaries +from bird.preProcess.stl_patch.stl_shapes import * + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser(description="Generate boundary patch") + parser.add_argument( + "-i", + "--input", + type=str, + metavar="", + required=False, + help="Boundary patch Json input", + default=os.path.join( + BIRD_PRE_PATCH_TEMP_DIR, "spider_spg/inlets_outlets.json" + ), + ) + parser.add_argument( + "-v", "--verbose", action="store_true", help="plot on screen" + ) + args = parser.parse_args() + bc_patch_dict = parseJsonFile(args.input) + write_boundaries(bc_patch_dict) + + if args.verbose: + # plot + from bird.utilities.stl_plotting import plotSTL, plt, pretty_labels + + axes = plotSTL("inlets.stl") + pretty_labels("x", "y", zlabel="z", fontsize=14, ax=axes) + plt.show() diff --git a/bird/__init__.py b/bird/__init__.py index dcff4926..84cc8487 100644 --- a/bird/__init__.py +++ b/bird/__init__.py @@ -7,6 +7,7 @@ BIRD_DIR = os.path.dirname(os.path.realpath(__file__)) BIRD_MESH_DIR = os.path.join(BIRD_DIR, "meshing") BIRD_POST_DIR = os.path.join(BIRD_DIR, "postProcess") +BIRD_PRE_DIR = os.path.join(BIRD_DIR, "preProcess") BIRD_BLOCK_CYL_MESH_TEMP_DIR = os.path.join( BIRD_MESH_DIR, "block_cyl_mesh_templates" ) @@ -25,5 +26,8 @@ BIRD_STIR_TANK_CASE_TEMP_DIR = os.path.join( BIRD_MESH_DIR, "stir_tank_case_templates" ) +BIRD_PRE_PATCH_TEMP_DIR = os.path.join( + BIRD_PRE_DIR, "stl_patch", "bc_patch_mesh_template" +) BIRD_EARLY_PRED_DATA_DIR = os.path.join(BIRD_POST_DIR, "data_early") BIRD_INV_DIR = os.path.join(BIRD_DIR, "inverse_modeling") diff --git a/bird/meshing/_stir_tank_reactor.py b/bird/meshing/_stir_tank_reactor.py index 3d988d47..b0dc13ae 100644 --- a/bird/meshing/_stir_tank_reactor.py +++ b/bird/meshing/_stir_tank_reactor.py @@ -6,7 +6,6 @@ class StirTankReactor: - def __init__( self, Dt, diff --git a/bird/meshing/block_rect_mesh.py b/bird/meshing/block_rect_mesh.py index d6a04a45..c594ddf6 100644 --- a/bird/meshing/block_rect_mesh.py +++ b/bird/meshing/block_rect_mesh.py @@ -42,7 +42,6 @@ def write_vertices(outfile, lengths, nboxes): def write_this_block(outfile, comment, ids, mesh, zonename="none"): - outfile.write("\n //" + comment + "\n") outfile.write("hex (") for i in range(len(ids)): @@ -57,7 +56,6 @@ def write_this_block(outfile, comment, ids, mesh, zonename="none"): def write_blocks(outfile, blockids, lengths, nboxes, points_per_len): - nx = nboxes[0] ny = nboxes[1] nz = nboxes[2] @@ -74,7 +72,6 @@ def write_blocks(outfile, blockids, lengths, nboxes, points_per_len): outfile.write("blocks\n(\n") for bl in range(nblocks): - i = blockids[bl][0] j = blockids[bl][1] k = blockids[bl][2] @@ -98,7 +95,6 @@ def write_blocks(outfile, blockids, lengths, nboxes, points_per_len): def write_patches(outfile): - outfile.write( "\n// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n" ) @@ -212,7 +208,6 @@ def assemble_mesh(input_file): def writeBlockMeshDict(out_folder, geomDict, meshDict): - n_blocks = len(geomDict["fluid_blocks"]) blockids = np.array(geomDict["fluid_blocks"]).reshape(n_blocks, 3) diff --git a/bird/meshing/stir_tank_mesh.py b/bird/meshing/stir_tank_mesh.py index daf89122..ac8535fd 100644 --- a/bird/meshing/stir_tank_mesh.py +++ b/bird/meshing/stir_tank_mesh.py @@ -11,7 +11,6 @@ def get_reactor_geom(yamlfile): def write_ofoam_preamble(outfile, react): - outfile.write( "/*--------------------------------*- C++ -*----------------------------------*\\\n" ) @@ -50,7 +49,6 @@ def write_ofoam_preamble(outfile, react): def write_vertices(outfile, react): - outfile.write( "\n// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n" ) @@ -67,7 +65,6 @@ def write_vertices(outfile, react): counter = 0 for repeat in range(2): for zi in range(nsections): - outfile.write("\n//section " + str(zi) + "\n") outfile.write("\n//center\n") outfile.write( @@ -120,7 +117,6 @@ def write_vertices(outfile, react): def get_globalindex_of(splti, ci, zi, react): - npts_per_section = react.npts_per_section centeroffset = react.centeroffset polyoffset = react.polyoffset @@ -140,7 +136,6 @@ def get_globalindex_of(splti, ci, zi, react): def get_baffle_point_of(splti, ci, zi, react): - baff_sections = react.baff_sections nsections = react.nsections npts_per_section = react.npts_per_section @@ -150,14 +145,11 @@ def get_baffle_point_of(splti, ci, zi, react): baffle_id = get_globalindex_of(splti, ci, zi, react) if zi in baff_sections: - if ci == hub_circ: - if splti % 2 == 0: # even number for impeller baffle_id += nsections * npts_per_section if ci == tank_circ: - if splti % 2 == 1: # odd number for baffle baffle_id += nsections * npts_per_section @@ -165,7 +157,6 @@ def get_baffle_point_of(splti, ci, zi, react): def write_edges(outfile, react): - nsections = react.nsections nsplits = react.nsplits dangle = react.dangle @@ -179,7 +170,6 @@ def write_edges(outfile, react): outfile.write("edges\n(\n") for zi in range(nsections): - outfile.write("\n//section " + str(zi) + "\n") offset = 1 + nsplits # one for center and nsplits for polygon @@ -212,7 +202,6 @@ def write_edges(outfile, react): def write_this_block(outfile, comment, ids, mesh, zonename="none"): - outfile.write("\n //" + comment + "\n") outfile.write("hex (") for i in range(len(ids)): @@ -227,7 +216,6 @@ def write_this_block(outfile, comment, ids, mesh, zonename="none"): def write_blocks(outfile, react): - outfile.write( "\n// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n" ) @@ -253,7 +241,6 @@ def write_blocks(outfile, react): hub_volumes = react.hub_volumes for zi in range(nvolumes): - outfile.write("\n//section " + str(zi) + "-" + str(zi + 1) + "\n") offset0 = zi * npts_per_section @@ -266,7 +253,6 @@ def write_blocks(outfile, react): # skip polygon blocks in stem sections if zi in nonstem_volumes: for i in range(nsplits): - localind1 = centeroffset + i % nsplits localind2 = centeroffset + (i + 1) % nsplits @@ -296,7 +282,6 @@ def write_blocks(outfile, react): outfile.write("\n//circles\n") for ci in range(ncirc): - zonename = "none" # skip blocks inside hub @@ -311,7 +296,6 @@ def write_blocks(outfile, react): zonename = "rotor" for i in range(nsplits): - idarray[0] = get_baffle_point_of(i, ci, zi, react) idarray[1] = get_baffle_point_of(i, ci, zi + 1, react) idarray[2] = get_globalindex_of(i + 1, ci, zi + 1, react) @@ -335,7 +319,6 @@ def write_blocks(outfile, react): def write_patches(outfile, react): - inhub_ci = react.inhub_circ hub_ci = react.hub_circ rot_ci = react.rot_circ @@ -454,7 +437,6 @@ def write_patches(outfile, react): zi_top = zi_bottom + 1 # bottom of impeller section for zi in [zi_bottom, zi_top]: - outfile.write("\n\t\t//hub to blade circle\n") for i in range(nsplits): outfile.write("\t\t( ") @@ -544,7 +526,6 @@ def write_patches(outfile, react): zi_pairs = [[hub_vol - 1, hub_vol], [hub_vol + 1, hub_vol + 2]] for zi_pair in zi_pairs: - zi_below = zi_pair[0] zi_above = zi_pair[1] @@ -685,7 +666,6 @@ def write_patches(outfile, react): zi_pairs.append([vols, vols + 1]) for zi_pair in zi_pairs: - zi_below = zi_pair[0] zi_above = zi_pair[1] outfile.write( @@ -712,7 +692,6 @@ def write_patches(outfile, react): outfile.write("\n\tempty inside_to_hub_copy\n\t(\n") for zi_pair in zi_pairs: - zi_below = zi_pair[0] zi_above = zi_pair[1] outfile.write( @@ -740,7 +719,6 @@ def write_patches(outfile, react): outfile.write("\n\tempty hub_to_rotor\n\t(\n") for zi_pair in zi_pairs: - zi_below = zi_pair[0] zi_above = zi_pair[1] outfile.write( @@ -766,7 +744,6 @@ def write_patches(outfile, react): outfile.write("\n\tempty hub_to_rotor_copy\n\t(\n") for zi_pair in zi_pairs: - zi_below = zi_pair[0] zi_above = zi_pair[1] outfile.write( diff --git a/bird/meshing/stl_mesh_tools.py b/bird/meshing/stl_mesh_tools.py deleted file mode 100644 index 1cf5546a..00000000 --- a/bird/meshing/stl_mesh_tools.py +++ /dev/null @@ -1,166 +0,0 @@ -import sys - -import numpy as np -import stl -from scipy.spatial import Delaunay - - -def tri_area(point1, point2, point3): - return 0.5 * ( - point1[0] * (point2[1] - point3[1]) - + point2[0] * (point3[1] - point1[1]) - + point3[0] * (point1[1] - point2[1]) - ) - - -def patch_area(tri): - vertices = tri.simplices - points = tri.points - area = 0 - for i_triangle in range(vertices.shape[0]): - area += tri_area( - points[vertices[i_triangle, 0], :], - points[vertices[i_triangle, 1], :], - points[vertices[i_triangle, 2], :], - ) - return area - - -def triangulate(vertices): - points = np.zeros((vertices.shape[0], 2)) - points[:, 0] = vertices[:, 0] - points[:, 1] = vertices[:, 2] - tri = Delaunay(points) - area = patch_area(tri) - return np.array(tri.simplices), area - - -def makePolygon(rad, nvert): - theta = 2 * np.pi / nvert - vertices = [] - for i in range(nvert): - vertices.append( - [ - rad * np.cos(theta * i + (np.pi / 2 - theta / 2)), - 0, - rad * np.sin(theta * i + (np.pi / 2 - theta / 2)), - ] - ) - vertices = np.array(vertices) - faces, area = triangulate(vertices) - - meshInpt = {} - meshInpt["vertices"] = vertices - meshInpt["faces"] = faces - - return meshInpt, area - - -def makeRectangle(w, h): - # Define vertices - vertices = np.array( - [ - [-w / 2, 0.0, h / 2], - [w / 2, 0.0, h / 2], - [w / 2, 0.0, -h / 2], - [-w / 2, 0.0, -h / 2], - ] - ) - - faces, area = triangulate(vertices) - - meshInpt = {} - meshInpt["vertices"] = vertices - meshInpt["faces"] = faces - - return meshInpt, area - - -def rotate(stlObj, theta=0): - stlObj.rotate([0, 1, 0], theta) - return stlObj - - -def translate(stlObj, vector=np.array([0, 0, 0])): - stlObj.translate(vector) - return stlObj - - -def traceMesh(meshInpt): - # Create the mesh - stlObj = stl.mesh.Mesh( - np.zeros(meshInpt["faces"].shape[0], dtype=stl.mesh.Mesh.dtype) - ) - for i, f in enumerate(meshInpt["faces"]): - for j in range(3): - stlObj.vectors[i][j] = meshInpt["vertices"][f[j], :] - - return stlObj - - -def makeSpider(centerRad, nArms, widthArms, lengthArms): - print("Making spider") - print(f"\tcenterRadius={centerRad}") - print(f"\tnArms={nArms}") - print(f"\twidthArms={widthArms}") - print(f"\tlengthArms={lengthArms}") - - globalArea = 0 - if nArms < 2: - print("ERROR: nArms must be >= 2") - print(f"Got nArms = {nArms}") - sys.exit() - if nArms == 2: - nVertPol = 4 - if nArms > 2: - nVertPol = nArms - centerMesh, centerArea = makePolygon(rad=centerRad, nvert=nVertPol) - globalArea += centerArea - vertices = centerMesh["vertices"] - maxWidth = np.linalg.norm((vertices[1, :] - vertices[0, :])) - if widthArms > maxWidth: - print("ERROR: arm width will make arms overlap") - print("Either increase center radius or reduce arm width") - sys.exit() - center = traceMesh(centerMesh) - - arms = [] - for i in range(nArms): - if nArms > 2: - if i < nArms - 1: - indp = i + 1 - indm = i - else: - indp = 0 - indm = i - if nArms == 2: - if i == 0: - indp = 1 - indm = 0 - else: - indp = 3 - indm = 2 - armMesh, armArea = makeRectangle(w=widthArms, h=lengthArms) - globalArea += armArea - arm = traceMesh(armMesh) - side = vertices[indp, :] - vertices[indm, :] - angle = np.arccos(np.dot(side, [1, 0, 0]) / np.linalg.norm(side)) - if side[2] <= 0: - angle *= -1 - arm = rotate(arm, angle) - midSide = (vertices[indp, :] + vertices[indm, :]) / 2 - trans = midSide * (1 + lengthArms / (2 * np.linalg.norm(midSide))) - arm = translate(arm, trans) - arms.append(arm) - - arms_data = [entry.data for entry in arms] - - combined = stl.mesh.Mesh(np.concatenate([center.data] + arms_data)) - - return combined, globalArea - - -def saveSTL(stlObj, filename="spg.stl"): - # Write the mesh to file - print(f"Generating {filename}") - stlObj.save(filename, mode=stl.Mode.ASCII) diff --git a/bird/postProcess/conditional_mean.py b/bird/postProcess/conditional_mean.py index f3db805c..19918882 100644 --- a/bird/postProcess/conditional_mean.py +++ b/bird/postProcess/conditional_mean.py @@ -17,7 +17,6 @@ def compute_cond_mean( diff_val_list=[], diff_name_list=[], ): - time_float_sorted, time_str_sorted = getCaseTimes(case_path) mesh_time_str = getMeshTime(case_path) cellCentres = readMesh( @@ -114,7 +113,6 @@ def sequencePlot( case_names=[], symbList=["-", "-d", "-^", "-.", "-s", "-o", "-+"], ): - if not len(case_names) == len(folder_names): case_names = [f"test{i}" for i in range(len(folder_names))] if len(case_names) > len(symbList): diff --git a/bird/preprocess/stl_patch/bc_patch_mesh_template/loop_reactor/inlets_outlets.json b/bird/preprocess/stl_patch/bc_patch_mesh_template/loop_reactor/inlets_outlets.json new file mode 100644 index 00000000..05018282 --- /dev/null +++ b/bird/preprocess/stl_patch/bc_patch_mesh_template/loop_reactor/inlets_outlets.json @@ -0,0 +1,11 @@ +{ + "inlets": [ + {"type": "circle", "centx": 5.0, "centy": 0.0, "centz": 0.5, "radius": 0.4, "normal_dir": 1,"nelements": 50}, + {"type": "circle", "centx": 2.5, "centy": 0.0, "centz": 0.5, "radius": 0.4, "normal_dir": 1,"nelements": 50}, + {"type": "circle", "centx": 7.5, "centy": 0.0, "centz": 0.5, "radius": 0.4, "normal_dir": 1,"nelements": 50} + ], + "outlets": [ + {"type": "circle", "centx": 0.5, "centy": 5.0, "centz": 0.5, "radius": 0.4, "normal_dir": 1,"nelements": 50}, + {"type": "circle", "centx": 0.5, "centy": 5.0, "centz": 4.5, "radius": 0.4, "normal_dir": 1,"nelements": 50} + ] +} diff --git a/bird/preprocess/stl_patch/bc_patch_mesh_template/spider_spg/inlets_outlets.json b/bird/preprocess/stl_patch/bc_patch_mesh_template/spider_spg/inlets_outlets.json new file mode 100644 index 00000000..16284d41 --- /dev/null +++ b/bird/preprocess/stl_patch/bc_patch_mesh_template/spider_spg/inlets_outlets.json @@ -0,0 +1,8 @@ +{ + "inlets": [ + {"type": "spider", "centerRad": 0.25, "nArms": 12, "widthArms": 0.1, "lengthArms": 0.5, "centx": 1, "centy": 0.0, "centz": 0, "normal_dir": 1} + ], + "outlets": [ + {"type": "circle", "centx": 0, "centy": 5.0, "centz": 0, "radius": 1, "normal_dir": 1,"nelements": 50} + ] +} diff --git a/bird/preprocess/stl_patch/stl_bc.py b/bird/preprocess/stl_patch/stl_bc.py new file mode 100644 index 00000000..4a51ce78 --- /dev/null +++ b/bird/preprocess/stl_patch/stl_bc.py @@ -0,0 +1,39 @@ +import json +import sys + +import numpy as np +import stl + +from bird.meshing._mesh_tools import parseJsonFile +from bird.preProcess.stl_patch.stl_shapes import * + + +def check_input(input_dict): + assert isinstance(input_dict, dict) + for bound in input_dict: + assert isinstance(input_dict[bound], list) + for patch in input_dict[bound]: + assert isinstance(patch, dict) + + +def get_all_vert_faces(input_dict, boundary_name): + patch_mesh_list = create_boundary_patch_list(input_dict, boundary_name) + boundary_mesh = STLMesh() + boundary_mesh.from_mesh_list(patch_mesh_list) + return boundary_mesh + + +def write_boundaries(input_dict): + check_input(input_dict) + for boundary_name in input_dict.keys(): + print(f"Making {boundary_name}") + boundary_mesh = get_all_vert_faces(input_dict, boundary_name) + print(f"\tArea {boundary_mesh.area} m2") + boundary_mesh.save(f"{boundary_name}.stl") + + +if __name__ == "__main__": + input_dict = parseJsonFile( + "bc_patch_mesh_template/loop_reactor/inlets_outlets.json" + ) + write_boundaries(input_dict) diff --git a/bird/preprocess/stl_patch/stl_mesh.py b/bird/preprocess/stl_patch/stl_mesh.py new file mode 100644 index 00000000..7a568c1c --- /dev/null +++ b/bird/preprocess/stl_patch/stl_mesh.py @@ -0,0 +1,151 @@ +import numpy as np +import stl +from scipy.spatial import Delaunay + + +class STLMesh: + def __init__( + self, + vertices=None, + faces=None, + area=None, + planar=None, + normal_dir=None, + ): + self.vertices = vertices + self.faces = faces + self.area = area + self.normal_dir = normal_dir + self.planar = planar + self.status = {} + self.update() + + def update_status(self): + if self.vertices is None: + self.status["vertices"] = False + elif self.vertices.shape[1] == 3 and len(self.vertices.shape) == 2: + self.status["vertices"] = True + else: + self.status["vertices"] = False + if self.faces is None: + self.status["faces"] = False + elif self.vertices.shape[1] == 3 and len(self.faces.shape) == 2: + self.status["faces"] = True + else: + self.status["faces"] = False + if self.area is None: + self.status["area"] = False + elif isinstance(self.area, float): + self.status["area"] = True + else: + self.status["area"] = False + if self.normal_dir is None: + self.status["normal_dir"] = False + elif isinstance(self.normal_dir, int): + self.status["normal_dir"] = True + else: + self.status["normal_dir"] = False + if self.planar is None: + self.status["planar"] = False + elif isinstance(self.planar, bool): + self.status["planar"] = True + else: + self.status["planar"] = False + + def update(self): + self.update_status() + if ( + self.status["vertices"] + and not self.status["normal_dir"] + and self.status["planar"] + ): + print("\t\tUpdating normal_dir to") + for i in range(3): + A = np.where(self.vertices[:, i] == self.vertices[0, i])[0] + if len(A) == self.vertices.shape[0]: + self.normal_dir = i + self.status["normal_dir"] = True + print(f"\t\t\t{i}") + break + if ( + self.status["vertices"] + and not self.status["faces"] + and self.status["planar"] + ): + print("\t\tUpdating faces") + points = np.zeros((self.vertices.shape[0], 2)) + count = 0 + for i in range(3): + if not i == self.normal_dir: + points[:, count] = self.vertices[:, i] + count += 1 + tri = Delaunay(points) + self.faces = np.array(tri.simplices) + self.status["faces"] = True + + if self.status["vertices"] and self.status["faces"]: + print("\t\tUpdating area") + self.calc_area() + self.status["area"] = True + + def to_stl(self): + stlObj = stl.mesh.Mesh( + np.zeros(self.faces.shape[0], dtype=stl.mesh.Mesh.dtype) + ) + for i, f in enumerate(self.faces): + for j in range(3): + stlObj.vectors[i][j] = self.vertices[f[j], :] + return stlObj + + def from_stl(self, stlObj): + self.vertices = np.unique(np.reshape(stlObj.vectors, (-1, 3)), axis=0) + self.faces = np.zeros((stlObj.vectors.shape[0], 3), dtype=int) + for i in range(stlObj.vectors.shape[0]): + for j in range(stlObj.vectors.shape[1]): + ind = np.argwhere( + np.linalg.norm( + self.vertices - stlObj.vectors[i, j], axis=1 + ) + == 0 + )[0][0] + self.faces[i, j] = ind + + def calc_area(self): + stlObj = self.to_stl() + stlObj.update_areas() + self.area = np.sum(stlObj.areas) + + def from_mesh_list(self, mesh_list): + if len(mesh_list) > 1: + self.planar = False + self.update_status() + + offset = int(0) + self.area = 0 + for patch_mesh in mesh_list: + patch_mesh.faces += offset + offset += len(patch_mesh.vertices) + if self.vertices is None: + self.vertices = patch_mesh.vertices + else: + self.vertices = np.vstack((self.vertices, patch_mesh.vertices)) + if self.faces is None: + self.faces = patch_mesh.faces + else: + self.faces = np.vstack((self.faces, patch_mesh.faces)) + self.area += patch_mesh.area + self.update() + + def rotate(self, theta=0, normal_dir=1): + stlObj = self.to_stl() + normal = [0, 0, 0] + normal[normal_dir] = 1 + stlObj.rotate(normal, theta) + return self.from_stl(stlObj) + + def translate(self, vector=np.array([0, 0, 0])): + self.vertices += vector + + def save(self, filename): + stlObj = self.to_stl() + stlObj.save(f"{filename}", mode=stl.Mode.ASCII) diff --git a/bird/preprocess/stl_patch/stl_shapes.py b/bird/preprocess/stl_patch/stl_shapes.py new file mode 100644 index 00000000..bfbaef1c --- /dev/null +++ b/bird/preprocess/stl_patch/stl_shapes.py @@ -0,0 +1,160 @@ +import json +import sys + +import numpy as np + +from bird.preProcess.stl_patch.stl_mesh import STLMesh + + +def make_polygon(rad, nvert, center, normal_dir): + print("\tMaking polygon") + theta = 2 * np.pi / nvert + vertices = np.zeros((nvert, 3)) + t1dir = (normal_dir + 1) % 3 + t2dir = (normal_dir + 2) % 3 + for i in range(nvert): + vertices[i, t2dir] = rad * np.cos(theta * i + (np.pi / 2 - theta / 2)) + vertices[i, t1dir] = rad * np.sin(theta * i + (np.pi / 2 - theta / 2)) + + stl_mesh = STLMesh(vertices=vertices, planar=True, normal_dir=normal_dir) + stl_mesh.translate(center) + + return stl_mesh + + +def make_rectangle(w, h, center, normal_dir): + print("\tMaking rectangle") + # Define vertices + t1dir = (normal_dir + 1) % 3 + t2dir = (normal_dir + 2) % 3 + vertices = np.zeros((4, 3)) + vertices[0, t2dir] = -w / 2 + vertices[0, t1dir] = h / 2 + vertices[1, t2dir] = w / 2 + vertices[1, t1dir] = h / 2 + vertices[2, t2dir] = w / 2 + vertices[2, t1dir] = -h / 2 + vertices[3, t2dir] = -w / 2 + vertices[3, t1dir] = -h / 2 + + stl_mesh = STLMesh(vertices=vertices, planar=True, normal_dir=normal_dir) + stl_mesh.translate(center) + + return stl_mesh + + +def make_circle(radius, center, normal_dir, npts=3): + print("\tMaking circle") + vertices = np.zeros((npts + 1, 3)) + t1dir = (normal_dir + 1) % 3 + t2dir = (normal_dir + 2) % 3 + for i in range(npts): + ang = i * 2.0 * np.pi / npts + vertices[i + 1, t2dir] = radius * np.cos(ang) + vertices[i + 1, t1dir] = radius * np.sin(ang) + + stl_mesh = STLMesh(vertices=vertices, planar=True, normal_dir=normal_dir) + stl_mesh.translate(center) + + return stl_mesh + + +def make_spider(centerRad, nArms, widthArms, lengthArms, center, normal_dir): + globalArea = 0 + if nArms < 2: + print("ERROR: nArms must be >= 2") + print(f"Got nArms = {nArms}") + sys.exit() + if nArms == 2: + nVertPol = 4 + if nArms > 2: + nVertPol = nArms + center_mesh = make_polygon( + rad=centerRad, nvert=nVertPol, center=(0, 0, 0), normal_dir=normal_dir + ) + vertices = center_mesh.vertices + maxWidth = np.linalg.norm((vertices[1, :] - vertices[0, :])) + if widthArms > maxWidth: + print("ERROR: arm width will make arms overlap") + print("Either increase center radius or reduce arm width") + sys.exit() + + arm_mesh_list = [] + for i in range(nArms): + if nArms > 2: + if i < nArms - 1: + indp = i + 1 + indm = i + else: + indp = 0 + indm = i + if nArms == 2: + if i == 0: + indp = 1 + indm = 0 + else: + indp = 3 + indm = 2 + arm_mesh = make_rectangle( + w=widthArms, h=lengthArms, center=(0, 0, 0), normal_dir=normal_dir + ) + side = vertices[indp, :] - vertices[indm, :] + if normal_dir == 1: + angle = np.arccos(np.dot(side, [1, 0, 0]) / np.linalg.norm(side)) + if side[2] <= 0: + angle *= -1 + else: + raise NotImplementedError + arm_mesh.rotate(angle, normal_dir=normal_dir) + midSide = (vertices[indp, :] + vertices[indm, :]) / 2 + trans = midSide * (1 + lengthArms / (2 * np.linalg.norm(midSide))) + arm_mesh.translate(trans) + arm_mesh_list.append(arm_mesh) + + spider_mesh = STLMesh() + spider_mesh.from_mesh_list([center_mesh] + arm_mesh_list) + spider_mesh.translate(center) + spider_mesh.planar = True + spider_mesh.normal_dir = normal_dir + + return spider_mesh + + +def create_boundary_patch_list(input_dict, boundary_name): + patch_mesh_list = [] + for patch in input_dict[boundary_name]: + if patch["type"].lower() == "circle": + patch_mesh = make_circle( + radius=patch["radius"], + center=(patch["centx"], patch["centy"], patch["centz"]), + normal_dir=patch["normal_dir"], + npts=patch["nelements"], + ) + elif patch["type"].lower() == "spider": + patch_mesh = make_spider( + centerRad=patch["centerRad"], + nArms=patch["nArms"], + widthArms=patch["widthArms"], + lengthArms=patch["lengthArms"], + center=(patch["centx"], patch["centy"], patch["centz"]), + normal_dir=patch["normal_dir"], + ) + elif patch["type"].lower() == "polygon": + patch_mesh = make_polygon( + rad=patch["radius"], + nvert=patch["npts"], + center=(patch["centx"], patch["centy"], patch["centz"]), + normal_dir=patch["normal_dir"], + ) + elif patch["type"].lower() == "rectangle": + patch_mesh = make_rectangle( + w=patch["width"], + h=patch["height"], + center=(patch["centx"], patch["centy"], patch["centz"]), + normal_dir=patch["normal_dir"], + ) + else: + raise NotImplementedError + patch_mesh_list.append(patch_mesh) + + return patch_mesh_list diff --git a/papers/co2_model/calibration/binBreakup/bcr_nn_list.py b/papers/co2_model/calibration/binBreakup/bcr_nn_list.py index 5ed35779..d01dc9ba 100644 --- a/papers/co2_model/calibration/binBreakup/bcr_nn_list.py +++ b/papers/co2_model/calibration/binBreakup/bcr_nn_list.py @@ -129,14 +129,19 @@ def make_data(self, np_file=None, scaler_folder="."): joblib.dump(scaler_y, os.path.join(scaler_folder, "scaler_y.mod")) # Split and shuffle - z_train, z_test, par_train, par_test, y_train, y_test = ( - train_test_split( - scaler_z.transform(z_dat), - scaler_par.transform(par_dat), - scaler_y.transform(y_dat), - test_size=0.01, - random_state=42, - ) + ( + z_train, + z_test, + par_train, + par_test, + y_train, + y_test, + ) = train_test_split( + scaler_z.transform(z_dat), + scaler_par.transform(par_dat), + scaler_y.transform(y_dat), + test_size=0.01, + random_state=42, ) return { diff --git a/papers/co2_model/calibration/breakup/bcr_nn_list.py b/papers/co2_model/calibration/breakup/bcr_nn_list.py index 5ed35779..d01dc9ba 100644 --- a/papers/co2_model/calibration/breakup/bcr_nn_list.py +++ b/papers/co2_model/calibration/breakup/bcr_nn_list.py @@ -129,14 +129,19 @@ def make_data(self, np_file=None, scaler_folder="."): joblib.dump(scaler_y, os.path.join(scaler_folder, "scaler_y.mod")) # Split and shuffle - z_train, z_test, par_train, par_test, y_train, y_test = ( - train_test_split( - scaler_z.transform(z_dat), - scaler_par.transform(par_dat), - scaler_y.transform(y_dat), - test_size=0.01, - random_state=42, - ) + ( + z_train, + z_test, + par_train, + par_test, + y_train, + y_test, + ) = train_test_split( + scaler_z.transform(z_dat), + scaler_par.transform(par_dat), + scaler_y.transform(y_dat), + test_size=0.01, + random_state=42, ) return { diff --git a/tests/meshing/test_stl.py b/tests/meshing/test_stl.py deleted file mode 100644 index 589c8de2..00000000 --- a/tests/meshing/test_stl.py +++ /dev/null @@ -1,22 +0,0 @@ -import sys - -import numpy as np - -from bird.meshing.stl_mesh_tools import makeSpider, saveSTL -from bird.utilities.stl_plotting import plotSTL, plt, pretty_labels - - -def test_spider(): - # Spider - combined, globalArea = makeSpider( - centerRad=0.25, - nArms=12, - widthArms=0.1, - lengthArms=0.5, - ) - print(f"\tglobalArea = {globalArea}") - saveSTL(combined, filename="spg.stl") - - # plot - axes = plotSTL("spg.stl") - pretty_labels("x", "y", zlabel="z", fontsize=14, ax=axes) diff --git a/tests/preProcess/test_stl_patch.py b/tests/preProcess/test_stl_patch.py new file mode 100644 index 00000000..ce17c479 --- /dev/null +++ b/tests/preProcess/test_stl_patch.py @@ -0,0 +1,29 @@ +import os + +import numpy as np + +from bird import BIRD_PRE_PATCH_TEMP_DIR +from bird.preProcess.stl_patch.stl_bc import parseJsonFile, write_boundaries +from bird.utilities.stl_plotting import plotSTL, pretty_labels + + +def test_spider_sparger(): + input_dict = parseJsonFile( + os.path.join(BIRD_PRE_PATCH_TEMP_DIR, "spider_spg/inlets_outlets.json") + ) + write_boundaries(input_dict) + # plot + axes = plotSTL("inlets.stl") + pretty_labels("x", "y", zlabel="z", fontsize=14, ax=axes) + + +def test_loop_reactor(): + input_dict = parseJsonFile( + os.path.join( + BIRD_PRE_PATCH_TEMP_DIR, "loop_reactor/inlets_outlets.json" + ) + ) + write_boundaries(input_dict) + # plot + axes = plotSTL("inlets.stl") + pretty_labels("x", "y", zlabel="z", fontsize=14, ax=axes)