Added LIFE project

joconnor22 · Oct 18, 2019 · 651732d · 651732d
commit 651732d
Show file tree

Hide file tree

Showing 37 changed files with 7,498 additions and 0 deletions.
diff --git a/LICENSE b/LICENSE
diff --git a/README.md b/README.md
@@ -0,0 +1,59 @@
+# LIFE
+### Overview
+A lattice Boltzmann-immersed boundary-finite element (LIFE) solver for fluid-structure interaction problems involving slender structures.
+
+### Key Features
+ - BGK lattice Boltzmann method (LBM)
+ - velocity and pressure boundary conditions set via the regularised method (Latt et al. 2008)
+ - convective outlet condition (Yang 2013)
+ - implicit immersed boundary method (IBM) for rigid and flexible structures (Li et al. 2016, Pinelli et al. 2010)
+ - flexible motion solved via a corotational finite element method (FEM) for slender structures
+ - implicit Newmark time integration scheme for FEM
+ - force/displacement mapping between LBM and FEM grids
+ - strongly-coupled block Gauss-Seidel implicit fluid-structure coupling scheme (Degroote 2013)
+ - dynamic under-relaxation with Aitken's method (Irons and Tuck 1969)
+ - shared memory parallelisation with OpenMP
+
+# Instruction for Users
+### Project Structure
+The LIFE project is organised into five folders:
+
+ - **examples**: contains the input files for some example cases
+ - **inc**: contains the project header files
+ - **input**: contains the geometry.config file that specifies the immersed boundary bodies
+ - **obj**: contains the intermediate object files that are generated during the build
+ - **src**: contains the project source files
+
+### Case Definition
+LIFE used two files to define the case setup:
+
+ - **params.h**: sets up the domain and flow conditions
+ - **geometry.config**: specifies the immersed boundary bodies (rigid/flexible)
+
+The **params.h** file is a header file and therefore the project must be rebuilt every time this file is modified. The **geometry.config** is read by LIFE during initialisation and therefore does not require a rebuild after every modification. To set up a case with no immersed boundary bodies then **geometry.config** can be left blank or removed entirely. If present, it must exist within the **input** folder within the working directory. The instructions for the case setup are given in the comments within each of these files.
+
+There are three steps to adding a new type of geometry for LIFE:
+
+ 1. Create a new keyword configuration describing the geometry in the **geometry.config** file.
+ 2. Add a new condition in the **ObjectsClass::geometryReadIn()** routine to read this configuration.
+ 3. Implement a new constructor for **IBMBodyClass** to build the new geometry type.
+
+### Dependencies
+LIFE requires LAPACK and the Boost library. To install these on Ubuntu type the following into a terminal:
+
+	sudo apt-get install libblas-dev liblapack-dev libboost-all-dev
+
+### Building and running LIFE
+A makefile is provided to build LIFE. To use the makefile just navigate to the top level LIFE directory and type:
+
+	make
+
+An executable named **LIFE** will be generated. To run this just type:
+
+	./LIFE
+
+### Post-Processing
+LIFE creates a directory called **Results**. This folder contains VTK and log files describing the results of the simulation.
+
+### Restarting
+LIFE has support for restarting a simulation from where it left off. To do this just re-run the executable. There is also support for changing the **geometry.config** between restarts to add/remove bodies mid-simulation.
diff --git a/examples/ChannelFlow/params.h b/examples/ChannelFlow/params.h
@@ -0,0 +1,111 @@
+/*
+    LIFE: Lattice boltzmann-Immersed boundary-Finite Element
+    Copyright (C) 2019 Joseph O'Connor
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef PARAMS_H	// PARAMS_H
+#define PARAMS_H
+
+// Include defs
+#include "defs.h"
+
+// Set number of OMP threads (if commented then it will use system max)
+//#define THREADS 12
+
+// Set resFactor (for easily changing mesh resolution)
+const int resFactor = 1;
+
+// Simulation options
+//#define INLET_RAMP 2.0			// Inlet velocity ramp-up
+//#define WOMERSLEY 5.0				// Womersley number for oscillating pressure gradients
+//#define UNI_EPSILON				// Calculate epsilon over all IBM bodies at once
+//#define ORDERED					// For deterministic reduction operations
+//#define INITIAL_DEFLECT 0.01		// Set an initial deflection (fraction of L)
+
+// Outputs
+#define VTK								// Write out VTK
+//#define VTK_FEM						// Write out the FEM VTK
+//#define FORCES						// Write out forces on structures
+//#define TIPS							// Write out tip positions
+
+// Domain setup (lattice)
+const int Nx = resFactor * 500 + 1;   	// Number of lattice sites in x-direction
+const int Ny = resFactor * 50 + 1;		// Number of lattice sites in y-direction
+
+// Domain setup (physical)
+const double height_p = 1.0;			// Domain height (m)
+const double rho_p = 1.0;				// Fluid density (kg/m^3)
+const double nu_p = 0.01;				// Fluid kinematic viscosity (m^2/s)
+
+// Initial conditions
+const double ux0_p = 0.0;				// Initial x-velocity (m/s)
+const double uy0_p = 0.0;				// Initial y-velocity (m/s)
+
+// Gravity and pressure gradient
+const double gravityX = 0.0;			// Gravity component in x-direction (m/s^2)
+const double gravityY = 0.0;			// Gravity component in y-direction (m/s^2)
+const double dpdx = 0.0;				// Pressure gradient in x-direction (Pa/m)
+const double dpdy = 0.0;				// Pressure gradient in x-direction (Pa/m)
+
+// Boundary conditions (set to eFluid for periodic)
+#define WALL_LEFT	eVelocity			// Boundary condition at left wall
+#define WALL_RIGHT	ePressure			// Boundary condition at right wall
+#define WALL_BOTTOM	eWall				// Boundary condition at bottom wall
+#define WALL_TOP	eWall				// Boundary condition at top wall
+//#define PROFILE 	eParabolic			// Profile shape (uncomment for uniform)
+
+// Inlet conditions
+const double uxInlet_p = 1.0;			// Inlet x-velocity (m/s)
+const double uyInlet_p = 0.0;			// Inlet y-velocity (m/s)
+
+// FEM Newmark integration parameters
+const double alpha = 0.25;				// Alpha parameter
+const double delta = 0.5;				// Delta parameter
+
+// FSI Coupling parameters
+const double relaxMax = 1.0;			// Relaxation parameter (initial guess)
+const double subTol = 1e-8;				// Tolerance for subiterations
+
+// Set omega by three different methods (comment/uncomment)
+// Set omega directly
+//const double omega = 1.875;
+//const double tStep = pow(1.0 / sqrt(3.0), 2.0) * pow(height_p / (Ny - 1), 2.0) * (1.0 - 0.5 * omega) / (nu_p * omega);
+
+// Set time step
+const double tStep = 0.0005 / (resFactor * resFactor);
+const double omega = 1.0 / (nu_p * tStep / (pow(1.0 / sqrt(3.0), 2.0) * pow(height_p / (Ny - 1), 2.0)) + 0.5);
+
+// Set lattice velocity
+//const double uLB = (2.0 / 3.0) * 0.2 * 1.0 / sqrt(3.0);
+//const double omega = 1.0 / (uLB * (Ny-1) / (uxInlet_p * height_p / (3.0 * nu_p)) + 0.5);
+
+// Number of time steps and how often to write out
+const int nSteps = static_cast<int>(round(30.0 / tStep));	// Number of timesteps
+const int tinfo = nSteps / 1000;							// Frequency to write out info and logs
+const int tVTK = nSteps / 100;								// Frequency to write out VTK
+const int tRestart = nSteps / 10;							// Frequency to write out restart files
+
+// Reference values
+const double ref_nu = nu_p;			// Reference kinematic viscosity (m^2/s)
+const double ref_rho = rho_p;		// Reference density (kg/m^3)
+const double ref_P = 0.0;			// Reference pressure (Pa)
+const double ref_L = height_p;		// Reference length (m)
+const double ref_U = uxInlet_p;		// Reference velocity (m/s)
+
+// Precision for ASCII output
+const int PRECISION = 10;
+
+#endif	// PARAMS_H
diff --git a/examples/Cylinder/input/geometry.config b/examples/Cylinder/input/geometry.config
@@ -0,0 +1,18 @@
+# This is the geometry configuration file. Each specific body is input via
+# each line within this file. This file should be placed within the input 
+# folder prior to running LIFE. The general format for specifying a body
+# is KEYWORD NUMBER ID POSITION SPACING OBJECT_SPECIFIC_PARAMETERS. An 
+# example for each keyword is given below:
+#
+#
+# CIRCLE:
+# CIRCLE NUMBER STARTID STARTCENTREX STARTCENTREY SPACING_X SPACING_Y RADIUS
+#
+#
+# FILAMENT:
+# FILAMENT NUMBER STARTID STARTX STARTY SPACING_X SPACING_Y LENGTH HEIGHT ANGLE FLEX_RIGID N_ELEMENTS BC DENSITY YOUNG_MODULUS
+#
+#
+#
+# ** START OF GEOMETRY CONFIG FILE: **
+CIRCLE 1 0 0.2 0.2 0.0 0.0 0.05
diff --git a/examples/Cylinder/params.h b/examples/Cylinder/params.h
@@ -0,0 +1,111 @@
+/*
+    LIFE: Lattice boltzmann-Immersed boundary-Finite Element
+    Copyright (C) 2019 Joseph O'Connor
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef PARAMS_H	// PARAMS_H
+#define PARAMS_H
+
+// Include defs
+#include "defs.h"
+
+// Set number of OMP threads (if commented then it will use system max)
+//#define THREADS 12
+
+// Set resFactor (for easily changing mesh resolution)
+const int resFactor = 2;
+
+// Simulation options
+#define INLET_RAMP 2.0				// Inlet velocity ramp-up
+//#define WOMERSLEY 5.0				// Womersley number for oscillating pressure gradients
+//#define UNI_EPSILON					// Calculate epsilon over all IBM bodies at once
+#define ORDERED						// For deterministic reduction operations
+//#define INITIAL_DEFLECT 0.01		// Set an initial deflection (fraction of L)
+
+// Outputs
+#define VTK								// Write out VTK
+//#define VTK_FEM						// Write out the FEM VTK
+#define FORCES							// Write out forces on structures
+//#define TIPS							// Write out tip positions
+
+// Domain setup (lattice)
+const int Nx = resFactor * 220 + 1;   	// Number of lattice sites in x-direction
+const int Ny = resFactor * 41 + 1;		// Number of lattice sites in y-direction
+
+// Domain setup (physical)
+const double height_p = 0.41;			// Domain height (m)
+const double rho_p = 1.0;				// Fluid density (kg/m^3)
+const double nu_p = 0.001;				// Fluid kinematic viscosity (m^2/s)
+
+// Initial conditions
+const double ux0_p = 0.0;				// Initial x-velocity (m/s)
+const double uy0_p = 0.0;				// Initial y-velocity (m/s)
+
+// Gravity and pressure gradient
+const double gravityX = 0.0;			// Gravity component in x-direction (m/s^2)
+const double gravityY = 0.0;			// Gravity component in y-direction (m/s^2)
+const double dpdx = 0.0;				// Pressure gradient in x-direction (Pa/m)
+const double dpdy = 0.0;				// Pressure gradient in x-direction (Pa/m)
+
+// Boundary conditions (set to eFluid for periodic)
+#define WALL_LEFT	eVelocity			// Boundary condition at left wall
+#define WALL_RIGHT	ePressure			// Boundary condition at right wall
+#define WALL_BOTTOM	eWall				// Boundary condition at bottom wall
+#define WALL_TOP	eWall				// Boundary condition at top wall
+#define PROFILE 	eParabolic			// Profile shape (uncomment for uniform)
+
+// Inlet conditions
+const double uxInlet_p = 1.0;			// Inlet x-velocity (m/s)
+const double uyInlet_p = 0.0;			// Inlet y-velocity (m/s)
+
+// FEM Newmark integration parameters
+const double alpha = 0.25;				// Alpha parameter
+const double delta = 0.5;				// Delta parameter
+
+// FSI Coupling parameters
+const double relaxMax = 1.0;			// Relaxation parameter (initial guess)
+const double subTol = 1e-8;				// Tolerance for subiterations
+
+// Set omega by three different methods (comment/uncomment)
+// Set omega directly
+//const double omega = 1.875;
+//const double tStep = pow(1.0 / sqrt(3.0), 2.0) * pow(height_p / (Ny - 1), 2.0) * (1.0 - 0.5 * omega) / (nu_p * omega);
+
+// Set time step
+const double tStep = 0.001 / (resFactor * resFactor);
+const double omega = 1.0 / (nu_p * tStep / (pow(1.0 / sqrt(3.0), 2.0) * pow(height_p / (Ny - 1), 2.0)) + 0.5);
+
+// Set lattice velocity
+//const double uLB = (2.0 / 3.0) * 0.2 * 1.0 / sqrt(3.0);
+//const double omega = 1.0 / (uLB * (Ny-1) / (uxInlet_p * height_p / (3.0 * nu_p)) + 0.5);
+
+// Number of time steps and how often to write out
+const int nSteps = static_cast<int>(round(20.0 / tStep));	// Number of timesteps
+const int tinfo = nSteps / 1000;							// Frequency to write out info and logs
+const int tVTK = nSteps / 200;								// Frequency to write out VTK
+const int tRestart = nSteps / 100;							// Frequency to write out restart files
+
+// Reference values
+const double ref_nu = nu_p;			// Reference kinematic viscosity (m^2/s)
+const double ref_rho = rho_p;		// Reference density (kg/m^3)
+const double ref_P = 0.0;			// Reference pressure (Pa)
+const double ref_L = 0.1;			// Reference length (m)
+const double ref_U = uxInlet_p;		// Reference velocity (m/s)
+
+// Precision for ASCII output
+const int PRECISION = 10;
+
+#endif	// PARAMS_H
diff --git a/examples/Honami/input/geometry.config b/examples/Honami/input/geometry.config
@@ -0,0 +1,18 @@
+# This is the geometry configuration file. Each specific body is input via
+# each line within this file. This file should be placed within the input 
+# folder prior to running LIFE. The general format for specifying a body
+# is KEYWORD NUMBER ID POSITION SPACING OBJECT_SPECIFIC_PARAMETERS. An 
+# example for each keyword is given below:
+#
+#
+# CIRCLE:
+# CIRCLE NUMBER STARTID STARTCENTREX STARTCENTREY SPACING_X SPACING_Y RADIUS
+#
+#
+# FILAMENT:
+# FILAMENT NUMBER STARTID STARTX STARTY SPACING_X SPACING_Y LENGTH HEIGHT ANGLE FLEX_RIGID N_ELEMENTS BC DENSITY YOUNG_MODULUS
+#
+#
+#
+# ** START OF GEOMETRY CONFIG FILE: **
+FILAMENT 128 0 10.0 0.0 0.5 0.0 1.0 0.02 90.0 FLEXIBLE 20 CLAMPED 50.000000 269536.024785