stir tank example

README.md

@@ -3,14 +3,40 @@
 
 ## Installation
 
-```
+```bash
 conda create  --name brd python=3.10
 conda activate brd
 pip install -e .
 ```
 
 ## Meshing
 
+### Generate Stir tank mesh
+
+```bash
+inp=brd/meshing/stir_tank_mesh_templates/base_tank/tank_par.yaml
+out=brd/meshing/stir_tank_case_templates/base/system/blockMeshDict
+
+python applications/write_stir_tank_mesh.py -i $inp -o $out
+```
+
+Generates a blockMeshDict
+
+Then activate openFoam environement (tested with OpenFoam9) and mesh with
+
+```bash
+blockMesh -dict system/blockMeshDict
+stitchMesh -perfect -overwrite inside_to_hub inside_to_hub_copy
+stitchMesh -perfect -overwrite hub_to_rotor hub_to_rotor_copy
+transformPoints "rotate=((0 0 1)(0 1 0))":
+```
+Mesh visualized in Paraview
+
+<p float="left">
+  <img src="assets/stir_tank.png" width="350"/>
+</p>
+
+
 ### Generate STL mesh
 
 `python applications/write_stl_mesh.py -v -cr 0.25 -na 12 -aw 0.1 -al 0.5`

applications/write_stir_tank_mesh.py

@@ -0,0 +1,39 @@
+import os
+import sys
+import numpy as np
+from brd.meshing.stir_tank_mesh import get_reactor_geom, write_ofoam_preamble, write_vertices, write_edges, write_blocks, write_patches
+from brd import BRD_STIR_TANK_MESH_TEMP_DIR
+import argparse
+
+def main():
+    parser = argparse.ArgumentParser(description="Stir tank meshing")
+    parser.add_argument(
+        "-o",
+        "--output_file",
+        type=str,
+        metavar="",
+        required=True,
+        help="Output blockMeshDict",
+        default="blockMeshDict",
+    )
+    parser.add_argument(
+        "-i",
+        "--react_in",
+        type=str,
+        metavar="",
+        required=True,
+        help="YAML file containing geometry details of reactor",
+        default=os.path.join(BRD_STIR_TANK_MESH_TEMP_DIR, "base_tank", "tank_par.yaml"),
+    )
+    args = parser.parse_args()
+    with open(args.output_file, "w") as outfile:
+        react = get_reactor_geom(args.react_in)
+        write_ofoam_preamble(outfile, react)
+        write_vertices(outfile, react)
+        write_edges(outfile, react)
+        write_blocks(outfile, react)
+        write_patches(outfile, react)
+
+
+if __name__ == "__main__":
+    main()

brd/__init__.py

@@ -6,4 +6,6 @@
 BRD_MESH_DIR = os.path.join(BRD_DIR, "meshing")
 BRD_BLOCK_CYL_MESH_TEMP_DIR = os.path.join(BRD_MESH_DIR, "block_cyl_mesh_templates")
 BRD_BLOCK_CYL_CASE_TEMP_DIR = os.path.join(BRD_MESH_DIR, "block_cyl_case_templates")
+BRD_STIR_TANK_MESH_TEMP_DIR = os.path.join(BRD_MESH_DIR, "stir_tank_mesh_templates")
+BRD_STIR_TANK_CASE_TEMP_DIR = os.path.join(BRD_MESH_DIR, "stir_tank_case_templates")
 BRD_INV_DIR = os.path.join(BRD_DIR, "inverse_modeling")

brd/meshing/_stir_tank_reactor.py

@@ -0,0 +1,110 @@
+from pathlib import Path
+from ruamel.yaml import YAML
+import os
+import numpy as np
+
+class StirTankReactor:
+
+    def __init__(self, Dt, Da, H, nimpellers, C, W, L, Lin, J, Wh, polyrad, Z0, nr, nz, Npoly, Na, nbaffles):
+        # Loop through params and setattr v to self.k
+        for k,v in locals().items():
+            if k!='self':
+                setattr(self, k, v)
+        self.Dh = Da - 2 * L
+        self.Dmrf = (Da+Dt-2*J)/2
+        self.nsplits=2*nbaffles    #we need twice the number of splits
+        self.dangle=2.0*np.pi/float(self.nsplits)
+
+        self.circradii=np.array([self.Dh/2-Lin,self.Dh/2,Da/2,self.Dmrf/2,Dt/2-J,Dt/2])
+        self.ncirc = len(self.circradii)
+        self.hub_circ   = 1
+        self.inhub_circ = self.hub_circ-1  #circle inside hub
+        self.rot_circ   = self.hub_circ+1
+        self.mrf_circ   = self.rot_circ+1
+        self.tank_circ = self.ncirc-1
+
+        self.reacthts = [Z0]
+        self.baff_sections = []
+        self.baff_volumes = []
+        self.hub_volumes=[]
+        count=1
+        for n_imp in range(self.nimpellers):
+            self.reacthts.append(Z0 + C[n_imp] -  W/2)
+
+            self.baff_sections.append(count)
+            self.baff_volumes.append(count)
+            count=count+1
+
+            self.reacthts.append(Z0 + C[n_imp] - Wh/2)
+
+            self.baff_sections.append(count)
+            self.baff_volumes.append(count)
+            self.hub_volumes.append(count)
+            count=count+1
+
+            self.reacthts.append(Z0 + C[n_imp] + Wh/2)
+
+            self.baff_sections.append(count)
+            self.baff_volumes.append(count)
+            count=count+1
+
+            self.reacthts.append(Z0 + C[n_imp] +  W/2)
+            self.baff_sections.append(count)
+            count=count+1
+
+        self.reacthts.append(Z0+H)
+
+
+        self.nsections = len(self.reacthts)
+        self.nvolumes = self.nsections-1
+        self.meshz = nz*np.diff(self.reacthts)
+        self.meshz = self.meshz.astype(int)+1 #avoid zero mesh elements
+
+        self.all_volumes=range(self.nvolumes)
+        self.nonbaff_volumes=[sec for sec in self.all_volumes if sec not in self.baff_volumes]
+        self.nonstem_volumes=[0,1] #this is 0,1 no matter how many impellers are there
+
+
+        #note: stem_volumes include hub volumes also
+        #these are volumes where we miss out polygon block
+        self.stem_volumes=[sec for sec in self.all_volumes if sec not in self.nonstem_volumes]
+
+        #removes hub_volumes here for declaring patches
+        self.only_stem_volumes=[sec for sec in self.stem_volumes if sec not in self.hub_volumes]
+
+        #to define mrf region
+        #not that [1] is not a stem volume but baffles are there
+        self.mrf_volumes=[1]+self.stem_volumes
+
+        #increase grid points in the impeller section
+        for i in self.baff_volumes:
+            self.meshz[i] *=2
+
+        self.meshr = nr*np.diff(self.circradii)
+
+        #adding polygon to hub mesh resolution
+        self.meshr = np.append(nr*polyrad,self.meshr)
+        self.meshr = self.meshr.astype(int)
+        self.meshr += 1 # to avoid being zero
+
+        self.centeroffset     = 1 #one point on the axis
+        self.polyoffset       = self.nsplits #number of points on polygon
+        self.npts_per_section = self.centeroffset + self.polyoffset + self.ncirc*self.nsplits #center+polygon+circles
+
+
+    @classmethod
+    def from_file(cls, yamlfile):
+        if '.yaml' not in yamlfile:
+            yamlfile += '.yaml'
+        if os.path.exists(yamlfile):
+            yamlpath = Path(yamlfile)
+        else:
+            raise FileNotFoundError(yamlfile)
+        yaml = YAML(typ='safe')
+        in_dict = yaml.load(yamlpath)
+        react_dict = {**in_dict["geometry"], **in_dict["mesh"]}
+        return cls(**react_dict)
+
+
+
+

brd/meshing/stir_tank_case_templates/base/constant/MRFProperties

@@ -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;
+    location    "constant";
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+MRF1
+{
+    cellZone    rotor;
+    active      yes;
+
+    // Fixed patches (by default they 'move' with the MRF zone)
+    nonRotatingPatches ();
+
+    origin    (0 0 0);
+    axis      (0 1 0);
+    omega     73.3; // rad/s  (1 rad/s = 9.55 rpm)
+}
+
+// ************************************************************************* //

brd/meshing/stir_tank_case_templates/base/constant/g

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  3.0.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    location    "constant";
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 -9.81 0);
+
+
+// ************************************************************************* //

brd/meshing/stir_tank_case_templates/base/constant/globalTransport

@@ -0,0 +1,37 @@
+//****** Temperature ************
+Temp   295.15; // K
+//*****Gas transport*******************
+muMixGas        1.479e-05;
+PrMixGas        0.7;
+LeGas     		1.0;
+dbubGas			0.001;
+//****Liquid transport*****************
+CpMixLiq        2359.8;
+muMixLiq        0.0029; // 1 cP at 20C, 0.467 at 60C 
+kThermLiq	0.167; // W/m-K
+RMixLiq         3000; //NOT SURE
+rho0MixLiq      810; // kg/m^3
+sigmaLiq        0.023; //surface tension N/m
+ddropLiq        0.003; // NOT SURE
+//Wilke-Chang params for diffusivity of aq. O2
+// not currently used, JJS 4/14/16
+WC_psi          2.6;
+WC_M            18; // kg/kmol
+WC_V            25.6e-3; // m3/kmol / O2=25.6e-3 / H2=14.3e-3 / CO2=34e-3
+//********************************************
+H_O2_298 0.032;
+DH_O2	1700;
+H_CO2_298 0.83;
+DH_CO2	2400;
+H_H2_298 0.019;
+DH_H2   500;
+H_CH4_298 0.032;
+DH_CH4   1900;
+// ******* Inlet conditions *********
+targetUs        0.02;
+inletArea       0.5;
+totalArea       1.0;
+inletAlphaGas   1.0;
+liqHeight       1.5;
+presfactor      1.1;  // (P1+rho*g*h)/(P1+rho*g*h/2)
+O2Liquid        0; // massFrac=C(mol/m3)*MW/rho 3.2e-6=0.1(mol/m3)