diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Allclean b/FSI/rigid_flap/OpenFOAM-RigidBody/Allclean new file mode 100755 index 000000000..be3a2d6f5 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Allclean @@ -0,0 +1,56 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory + +echo "Cleaning..." + +# Source tutorial clean functions +. $WM_PROJECT_DIR/bin/tools/CleanFunctions + +# Participant 1: Fluid (OpenFOAM) +Participant1="Fluid" +cd ${Participant1} + # Clean the case + cleanCase + rm -rfv 0 + # Create an empty .foam file for ParaView + # Note: ".foam" triggers the native OpenFOAM reader of ParaView. + # Change to ".OpenFOAM" to use the OpenFOAM reader provided with OpenFOAM. + touch ${Participant1}.foam +cd .. +# Remove the log files +rm -fv ${Participant1}_blockMesh.log +rm -fv ${Participant1}_checkMesh.log +rm -fv ${Participant1}_decomposePar.log +rm -fv ${Participant1}.log +rm -fv ${Participant1}_reconstructPar.log + +# Participant 2: Solid (rigid body) +Participant2="Solid" +cd ./Solid/ + # Clean the case + echo "Cleaning Solid case" + rm -fv solution-*.vtk +cd .. + +# Remove the log files +rm -fv ${Participant2}.log + +# Remove the preCICE-related log files +echo "Deleting the preCICE log files..." +rm -fv \ + precice-*.log \ + precice-*-events.json + +# Output files for preCICE versions before 1.2: +rm -fv \ + iterations-${Participant1}.txt iterations-${Participant2}.txt \ + convergence-${Participant1}.txt convergence-${Participant2}.txt \ + Events-${Participant1}.log Events-${Participant2}.log \ + EventTimings-${Participant1}.log EventTimings-${Participant2}.log + +rm -rfv precice-run +rm -fv .${Participant1}-${Participant2}.address +rm -rf coupling-meshes + +echo "Cleaning complete!" +#------------------------------------------------------------------------------ diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/0.orig/U b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/0.orig/U new file mode 100644 index 000000000..69dfcf932 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/0.orig/U @@ -0,0 +1,54 @@ +FoamFile +{ + version 2.0; + format ascii; + class volVectorField; + object U; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 -1 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + + flap + { + type movingWallVelocity; + value uniform (0 0 0); + } + + top + { + type noSlip; + } + + bottom + { + type noSlip; + } + + inlet + { + type fixedValue; + value uniform (1 0 0); + } + + outlet + { + type zeroGradient; + } + + front + { + type empty; + } + back + { + type empty; + } +} + +// ************************************************************************* // diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/0.orig/p b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/0.orig/p new file mode 100644 index 000000000..90e9358f4 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/0.orig/p @@ -0,0 +1,52 @@ +FoamFile +{ + version 2.0; + format ascii; + class volScalarField; + object p; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 2 -2 0 0 0 0]; + +internalField uniform 0; + +boundaryField +{ + flap + { + type zeroGradient; + } + + top + { + type zeroGradient; + } + + bottom + { + type zeroGradient; + } + + inlet + { + type zeroGradient; + } + + outlet + { + type fixedValue; + value uniform 0; + } + + front + { + type empty; + } + back + { + type empty; + } +} + +// ************************************************************************* // diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/0.orig/pointDisplacement b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/0.orig/pointDisplacement new file mode 100644 index 000000000..e12eaa933 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/0.orig/pointDisplacement @@ -0,0 +1,53 @@ +FoamFile +{ + version 2.0; + format ascii; + class pointVectorField; + object pointDisplacement; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dimensions [0 1 0 0 0 0 0]; + +internalField uniform (0 0 0); + +boundaryField +{ + inlet + { + type fixedValue; + value uniform (0 0 0); + } + + outlet + { + type fixedValue; + value uniform (0 0 0); + } + + flap + { + type fixedValue; + value $internalField; + } + + top + { + type slip; + } + + bottom + { + type slip; + } + front + { + type empty; + } + back + { + type empty; + } +} + +// ************************************************************************* // diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/Fluid.foam b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/Fluid.foam new file mode 100644 index 000000000..e69de29bb diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/constant/dynamicMeshDict b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/constant/dynamicMeshDict new file mode 100644 index 000000000..d323bc059 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/constant/dynamicMeshDict @@ -0,0 +1,22 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object dynamicMeshDict; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +dynamicFvMesh dynamicMotionSolverFvMesh; + +motionSolverLibs ("libfvMotionSolvers.so"); + +solver displacementLaplacian; + +displacementLaplacianCoeffs { + + diffusivity quadratic inverseDistance (flap); + +} \ No newline at end of file diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/constant/transportProperties b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/constant/transportProperties new file mode 100644 index 000000000..870a0bb1c --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/constant/transportProperties @@ -0,0 +1,14 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location constant; + object transportProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + transportModel Newtonian; + + nu nu [0 2 -1 0 0 0 0 ] 0.001; + pRef pRef [1 -1 -2 0 0 0 0 ] 0.0; diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/constant/turbulenceProperties b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/constant/turbulenceProperties new file mode 100644 index 000000000..f906d4153 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/constant/turbulenceProperties @@ -0,0 +1,11 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location constant; + object turbulenceProperties; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +simulationType laminar; diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/blockMeshDict b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/blockMeshDict new file mode 100644 index 000000000..c8a955624 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/blockMeshDict @@ -0,0 +1,199 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +scale 1; + +b -0.1; // z-back +f 0.1; // z-front + +// Vertices +X1 -0.2; // pre rigid body +X2 0.0; // begin rigid body +X3 0.2; // end rigid body +X4 0.6; // wake + +Y1 -0.13; // distance to bottom +Y2 -0.01; // half body height +Y3 0.01; // half body height +Y4 0.12; // distance to top + +// Blocks +H1 12; +H2 14; +H3 20; + +V1 12; +V2 2; +V3 11; + +G3 4; + +vertices +( + // X1 layer back + ($X1 $Y1 $b) // 0 + ($X1 $Y2 $b) + ($X1 $Y3 $b) + ($X1 $Y4 $b) + + // X1 layer front + ($X1 $Y1 $f) // 4 + ($X1 $Y2 $f) + ($X1 $Y3 $f) + ($X1 $Y4 $f) + + // X2 layer back + ($X2 $Y1 $b) // 8 + ($X2 $Y2 $b) + ($X2 $Y3 $b) + ($X2 $Y4 $b) + + // X2 layer front + ($X2 $Y1 $f) // 12 + ($X2 $Y2 $f) + ($X2 $Y3 $f) + ($X2 $Y4 $f) + + // X3 layer back + ($X3 $Y1 $b) // 16 + ($X3 $Y2 $b) + ($X3 $Y3 $b) + ($X3 $Y4 $b) + + // X3 layer front + ($X3 $Y1 $f) // 20 + ($X3 $Y2 $f) + ($X3 $Y3 $f) + ($X3 $Y4 $f) + + // X4 layer back + ($X4 $Y1 $b) // 24 + ($X4 $Y2 $b) + ($X4 $Y3 $b) + ($X4 $Y4 $b) + + // X4 layer front + ($X4 $Y1 $f) // 28 + ($X4 $Y2 $f) + ($X4 $Y3 $f) + ($X4 $Y4 $f) + +); + +blocks +( + // Block 0-3 + hex ( 0 8 9 1 4 12 13 5) ($H1 $V1 1) simpleGrading (1 1 1) + hex ( 1 9 10 2 5 13 14 6) ($H1 $V2 1) simpleGrading (1 1 1) + hex ( 2 10 11 3 6 14 15 7) ($H1 $V3 1) simpleGrading (1 1 1) + + // Block 4-6 \5 + hex ( 8 16 17 9 12 20 21 13) ($H2 $V1 1) simpleGrading (1 1 1) + hex ( 10 18 19 11 14 22 23 15) ($H2 $V3 1) simpleGrading (1 1 1) + + // Block 7-9 + hex ( 16 24 25 17 20 28 29 21) ($H3 $V1 1) simpleGrading ($G3 1 1) + hex ( 17 25 26 18 21 29 30 22) ($H3 $V2 1) simpleGrading ($G3 1 1) + hex ( 18 26 27 19 22 30 31 23) ($H3 $V3 1) simpleGrading ($G3 1 1) + +); + + +boundary +( + back + { + type empty; + faces + ( + (0 8 9 1) + (1 9 10 2) + (2 10 11 3) + (8 16 17 9) + (10 18 19 11) + (16 24 25 17) + (17 25 26 18) + (18 26 27 19) + ); + } + + front + { + type empty; + faces + ( + (4 12 13 5) + (5 13 14 6) + (6 14 15 7) + (12 20 21 13) + (14 22 23 15) + (20 28 29 21) + (21 29 30 22) + (22 30 31 23) + ); + } + + inlet + { + type patch; + faces + ( + (0 4 5 1) + (1 5 6 2) + (2 6 7 3) + ); + } + + outlet + { + type patch; + faces + ( + (24 28 29 25) + (25 29 30 26) + (26 30 31 27) + ); + } + + flap + { + type wall; + faces + ( + (9 13 14 10) + (9 17 21 13) + (10 18 22 14) + (18 22 21 17) + ); + } + + bottom + { + type wall; + faces + ( + (0 8 12 4) + (8 16 20 12) + (16 24 28 20) + ); + } + + top + { + type wall; + faces + ( + (3 11 15 7) + (11 19 23 15) + (19 27 31 23) + ); + } +); + +// ************************************************************************* // diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/controlDict b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/controlDict new file mode 100644 index 000000000..8a241d19b --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/controlDict @@ -0,0 +1,61 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location system; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// Note: With OpenFOAM v1806 and OpenFOAM 6, the DyM solvers +// were marked deprecated and merged into their respective standard solvers. +application pimpleFoam; // OpenFOAM v1806, OpenFOAM 6, or newer +// application pimpleDyMFoam; // OpenFOAM v1712, OpenFOAM 5.x, or older + + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 2.5; + +deltaT 2.5e-2; + +writeControl adjustableRunTime; + +writeInterval 2.5e-2; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 10; + +writeCompression off; + +timeFormat general; + +timePrecision 8; + +functions +{ + forces + { + type forces; + libs ( "libforces.so" ); + patches (flap); + rho rhoInf; + log true; + rhoInf 10; + CofR (0 0 0); + } + + preCICE_Adapter + { + type preciceAdapterFunctionObject; + libs ("libpreciceAdapterFunctionObject.so"); + } +} diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/decomposeParDict b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/decomposeParDict new file mode 100644 index 000000000..dee82bfd0 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/decomposeParDict @@ -0,0 +1,25 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location system; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + numberOfSubdomains 2; + + method hierarchical; + hierarchicalCoeffs + { + n (2 1 1); + delta 0.001; + order xyz; + } + + distributed false; + roots + ( + ); + diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/fvSchemes b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/fvSchemes new file mode 100644 index 000000000..c784ff608 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/fvSchemes @@ -0,0 +1,41 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + ddtSchemes + { + default backward; + } + + gradSchemes + { + default cellLimited Gauss linear 1; + } + + divSchemes + { + default none; + div(phi,U) Gauss linearUpwind grad(U); + div((nuEff*dev2(T(grad(U))))) Gauss linear; + } + + interpolationSchemes + { + default linear; + } + + laplacianSchemes + { + default Gauss linear corrected; + } + + snGradSchemes + { + default corrected; + } diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/fvSolution b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/fvSolution new file mode 100644 index 000000000..5592bdf82 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/fvSolution @@ -0,0 +1,85 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +solvers +{ + p + { + solver GAMG; + tolerance 1e-6; + relTol 1e-4; + smoother DICGaussSeidel; + } + + pFinal + { + $p; + tolerance 1e-07; + relTol 0; + } + + pcorr + { + solver GAMG; + tolerance 1e-5; + relTol 1e-3; + smoother GaussSeidel; + } + + pcorrFinal + { + $pcorr; + relTol 0; + } + + phi + { + $p; + } + + "(U|cellDisplacement)" + { + solver smoothSolver; + smoother symGaussSeidel; + tolerance 1e-7; + relTol 1e-5; + } + + "(U|cellDisplacement)Final" + { + $U; + relTol 0; + } +} + +PIMPLE +{ + nOuterCorrectors 10; + nCorrectors 2; + nNonOrthogonalCorrectors 1; + tolerance 1.0e-8; + + correctPhi yes; + relTol 1e-4; + pisoTol 1e-6; + consistent true; + +} +PISO +{ + nNonOrthogonalCorrectors 1; +} +potentialFlow +{ + nNonOrthogonalCorrectors 1; +} + + +// ************************************************************************* // diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/preciceDict b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/preciceDict new file mode 100644 index 000000000..1f44d0e70 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Fluid/system/preciceDict @@ -0,0 +1,54 @@ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object preciceDict; +} + +preciceConfig "precice-config.xml"; + +participant Fluid; + +modules (FSI); + +interfaces +{ + Interface1 + { + mesh Fluid-Mesh-Centers; + patches (flap); + locations faceCenters; + + readData + ( + ); + + writeData + ( + Forces + ); + }; + + Interface2 + { + mesh Fluid-Mesh-Nodes; + patches (flap); + locations faceNodes; + + readData + ( + Displacements + ); + + writeData + ( + ); + }; +}; + +FSI +{ + rho rho [1 -3 0 0 0 0 0] 1000; +} diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/README.md b/FSI/rigid_flap/OpenFOAM-RigidBody/README.md new file mode 100644 index 000000000..019cd1b4e --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/README.md @@ -0,0 +1,28 @@ +# FSI tutorial of a rigid body flap in a channel flow + +This tutorial deals with a fluid-structure-interaction problem. The fluid part of the simulation is computed using OpenFOAM and the rigid body motion is a rigid body model written in `c++` with only a single degree of freedom, namely the deflection angle of the flap in the channel. The rigid body is fixed in the origin at (0,0) and the force exerted by the fluid on the rigid body structure causes an oscillatory rotation of the body. The simulation runs for 2.5 seconds. After 1.5 seconds we increase the angular stiffness of the structure by a factor of 8 to stabilize the coupled problem. Feel free to modify these parameters and increase the simulation time. + +![overview](overview.png) + +## Building the `rigid_body` participant + +Before starting the coupled simulation, the `rigid_body` solver needs to be compiled using `cmake`. You can run the following commands from this directory to build the `rigid_body` solver +``` +cd Solid && cmake . && make +``` + +## Running the coupled simulation + +You may run the two simulations in two different terminals and watch their output on the screen by using the scripts `runFluid` (or `runFluid -parallel`) and `runSolid`. You can cleanup the simulation using `Allclean`. + +There is an [open issue](https://github.com/precice/openfoam-adapter/issues/26) that leads to additional "empty" result directories when running with some OpenFOAM versions, leading to inconveniences during post-processing. Please run the script `removeObsoleteSolvers.sh` to delete the additional files. + +Ini serial, the simulation takes roughly 30 seconds to compute. + +## Visualizing the results + +You can visualize the simulation results of the `Fluid` participant using paraView (use paraFoam to trigger the OpenFOAM native reader or load the (empty) file `Fluid.foam` into paraView). The rigid body doesn't generate any readable output files, but the motion can be observed in the OpenFOAM data. In addition, one could visualize the coupling meshes including the exchanged coupling data. preCICE generates the relevant files during the simulation and stores them in the directory `coupling-meshes`. In order to visualize the results, load the `vtk` files in paraView and apply a `Glyph` filter. Depending on the specific paraView version, you might additionally need to disable the `ScaleArray` option by selecting `No scale array` since the exchanged data might be inappropriate for a scaling operation. You can further add a `Warp By Vector` filter with `Displacements` to deform the coupling data. + +## Disclaimer + +This offering is not approved or endorsed by OpenCFD Limited, producer and distributor of the OpenFOAM software via www.openfoam.com, and owner of the OPENFOAM® and OpenCFD® trade marks. diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Solid/CMakeLists.txt b/FSI/rigid_flap/OpenFOAM-RigidBody/Solid/CMakeLists.txt new file mode 100644 index 000000000..4024b0f9d --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Solid/CMakeLists.txt @@ -0,0 +1,10 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 3.10.2) +SET(TARGET "rigid_body_solver") +PROJECT(${TARGET} LANGUAGES CXX DESCRIPTION "rigid_body_solver") + +FIND_PACKAGE(precice REQUIRED CONFIG) +ADD_EXECUTABLE( + ${TARGET} + ${TARGET}.cpp) + +TARGET_LINK_LIBRARIES(${TARGET} PRIVATE precice::precice) diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/Solid/rigid_body_solver.cpp b/FSI/rigid_flap/OpenFOAM-RigidBody/Solid/rigid_body_solver.cpp new file mode 100644 index 000000000..60c23d7b6 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/Solid/rigid_body_solver.cpp @@ -0,0 +1,256 @@ +#include +#include +#include +#include + +using Vector = std::vector; + +struct DataContainer { + void save_old_state(const Vector &vertices, + const double &theta, + const double &theta_dot) + { + old_vertices = vertices; + old_theta = theta; + old_theta_dot = theta_dot; + } + + void reload_old_state(Vector &vertices, + double &theta, + double &theta_dot) const + { + vertices = old_vertices; + theta = old_theta; + theta_dot = old_theta_dot; + } + + Vector old_vertices; + double old_theta; + double old_theta_dot; +}; + +class Solver { +public: + Solver(const double moment_of_inertia) + : moment_of_inertia(moment_of_inertia) + { + } + + void + solve(const Vector &forces, + const Vector &initial_vertices, + Vector & vertices, + double & theta, + double & theta_dot, + const double spring_constant, + const double delta_t) const + { + // Compute total moment M = x^{n} x f^{n+1} + double moment = 0; + for (unsigned int i = 0; i < forces.size() / 2; ++i) + moment += vertices[2 * i] * forces[2 * i + 1] - vertices[2 * i + 1] * forces[2 * i]; + + // Store rigid body angle at the previous time level theta^{n} + const double theta_old = theta; + + // Update angle to theta^{n+1} according to forward Euler method (simplified moment + // computation, which does not depend on the updated configuration). The contribution + // of the spring with strength k is discretized implicitly. + // theta^{n+1} = 1/ (1 - (k/I)*dt^2) * dt^2 * M / I + dt * \dot{theta}^{n} + theta^{n} + theta = (1. / (1 - (spring_constant / moment_of_inertia) * std::pow(delta_t, 2))) * + (std::pow(delta_t, 2) * moment / moment_of_inertia + delta_t * theta_dot + theta); + + // Update angular velocity + // \dot{theta}^{n+1} = (theta^{n+1} - \dot{theta}^{n}) / dt + theta_dot = (theta - theta_old) / delta_t; + + // Update vertices according to rigid body rotation using an out-of-plane (z-axis) rotation matrix + // x^{n+1} = x^{0} * cos(theta^{n+1}) + y^{0} * sin(theta^{n+1}) + // y^{n+1} = -x^{0} * sin(theta^{n+1}) + y^{0} * cos(theta^{n+1}) + for (uint i = 0; i < vertices.size() / 2; ++i) { + const double x_coord = initial_vertices[2 * i]; + vertices[2 * i] = x_coord * std::cos(theta) + initial_vertices[2 * i + 1] * std::sin(theta); + vertices[2 * i + 1] = -x_coord * std::sin(theta) + initial_vertices[2 * i + 1] * std::cos(theta); + } + std::cout << "Theta: " << theta << " Theta dot: " << theta_dot << " Moment: " << moment << " Spring force: " << spring_constant * theta << std::endl; + } + +private: + const double moment_of_inertia; +}; + +int main() +{ + std::cout << "Rigid body: starting... \n"; + + // Configuration settings + const std::string config_file_name("precice-config.xml"); + const std::string solver_name("Solid"); + const std::string mesh_name("Solid-Mesh"); + const std::string data_write_name("Displacements"); + const std::string data_read_name("Forces"); + + // Mesh configuration + constexpr int vertical_refinement = 3; + constexpr int horizontal_refinement = 6; + // Rotation centre is at (0,0) + constexpr double length = 0.2; + constexpr double height = 0.02; + + constexpr double density = 10000; + constexpr double spring_constant = -25; + + // Time, where spring is stiffened + constexpr double switch_time = 1.5; + constexpr double stiffening_factor = 8; + //*******************************************************************************************// + + // Derived quantities + constexpr int n_vertical_nodes = vertical_refinement * 2 + 1; + constexpr int n_horizontal_nodes = horizontal_refinement * 2 + 1; + // Substract shared nodes at each rigid body corner + constexpr int n_nodes = (n_vertical_nodes + n_horizontal_nodes - 2) * 2; + constexpr double mass = length * height * density; + // The moment of inertia is computed according to the rigid body configuration: + // a thin rectangular plate of height h, length l and mass m with axis of rotation + // at the end of the plate: I = (1/12)*m*(4*l^2+h^2) + constexpr double inertia_moment = (1. / 12) * mass * (4 * std::pow(length, 2) + std::pow(height, 2)); + constexpr double delta_y = height / (n_vertical_nodes - 1); + constexpr double delta_x = length / (n_horizontal_nodes - 1); + + // Create Solverinterface + precice::SolverInterface precice(solver_name, + config_file_name, + /*comm_rank*/ 0, + /*comm_size*/ 1); + + const int mesh_id = precice.getMeshID(mesh_name); + const int dim = precice.getDimensions(); + const int write_id = precice.getDataID(data_write_name, mesh_id); + const int read_id = precice.getDataID(data_read_name, mesh_id); + + // Set up data structures + Vector forces(dim * n_nodes); + Vector vertices(dim * n_nodes); + Vector displacement(dim * n_nodes); + std::vector vertex_ids(n_nodes); + double theta_dot = 0.0; + double theta = 0.0; + + { + // Define a boundary mesh + std::cout << "Rigid body: defining mesh: "; + std::cout << n_nodes << " nodes " << std::endl; + + // upper y + // o---o---o---o---o---o---o + // |.......................| + // lower x o.......................o upper x + // |.......................| + // o---o---o---o---o---o---o + // lower y + + // x planes + for (int i = 0; i < n_vertical_nodes; ++i) { + // lower x plane + vertices[dim * i] = 0.0; // fixed x + vertices[dim * i + 1] = -height * 0.5 + delta_y * i; // increasing y + // upper x plane + vertices[(2 * n_vertical_nodes) + dim * i] = length; // fixed x + vertices[(2 * n_vertical_nodes) + dim * i + 1] = vertices[dim * i + 1]; // increasing y + } + + // y planes + const unsigned int of = 2 * dim * n_vertical_nodes; // static offset + // Lower and upper bounds are already included due to positive/negative x-planes + const unsigned int n_remaining_nodes = n_horizontal_nodes - 2; + for (unsigned int i = 0; i < n_remaining_nodes; ++i) { + // lower y plane + vertices[of + dim * i] = delta_x * (i + 1); // increasing x + vertices[of + dim * i + 1] = -height * 0.5; // fixed y + // upper y plane + vertices[of + (2 * n_remaining_nodes) + dim * i] = vertices[of + dim * i]; // increasing x + vertices[of + (2 * n_remaining_nodes) + dim * i + 1] = -vertices[of + dim * i + 1]; // fixed y + } + } + // Store the initial configuration + const Vector initial_vertices = vertices; + + // Pass the vertices to preCICE + precice.setMeshVertices(mesh_id, + n_nodes, + vertices.data(), + vertex_ids.data()); + + // initialize the Solverinterface + double dt = precice.initialize(); + + // Compute the absolute displacement between the current vertices and the + // initial configuration. Here, this is mostly done for consistency reasons. + for (uint i = 0; i < displacement.size(); ++i) + displacement[i] = vertices[i] - initial_vertices[i]; + + // Set up a struct in order to store time dependent values + DataContainer data_container; + // Set up an object which handles the time integration + const Solver solver(inertia_moment); + + // Start time loop + double time = 0; + while (precice.isCouplingOngoing()) { + + std::cout << "Rigid body: t = " << time << "s \n"; + + std::cout << "Rigid body: reading initial data \n"; + if (precice.isReadDataAvailable()) + precice.readBlockVectorData(read_id, + n_nodes, + vertex_ids.data(), + forces.data()); + + // Store time dependent values + if (precice.isActionRequired( + precice::constants::actionWriteIterationCheckpoint())) { + data_container.save_old_state(vertices, theta, theta_dot); + + precice.markActionFulfilled( + precice::constants::actionWriteIterationCheckpoint()); + } + + const double current_spring = time > switch_time ? spring_constant * stiffening_factor : spring_constant; + // Solve system + solver.solve(forces, initial_vertices, vertices, theta, theta_dot, current_spring, dt); + + // Advance coupled system + // Compute absolute displacement with respect to the initial configuration + for (uint i = 0; i < displacement.size(); ++i) + displacement[i] = vertices[i] - initial_vertices[i]; + + std::cout << "Rigid body: writing coupling data \n"; + if (precice.isWriteDataRequired(dt)) + precice.writeBlockVectorData(write_id, + n_nodes, + vertex_ids.data(), + displacement.data()); + + std::cout << "Rigid body: advancing in time\n"; + dt = precice.advance(dt); + + // Reload time dependent values + if (precice.isActionRequired( + precice::constants::actionReadIterationCheckpoint())) { + data_container.reload_old_state(vertices, theta, theta_dot); + + precice.markActionFulfilled( + precice::constants::actionReadIterationCheckpoint()); + } + + // Increment time in case the time window has been completed + if (precice.isTimeWindowComplete()) + time += dt; + } + + std::cout << "Rigid body: closing...\n"; + + return 0; +} diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/overview.png b/FSI/rigid_flap/OpenFOAM-RigidBody/overview.png new file mode 100644 index 000000000..cf304937d Binary files /dev/null and b/FSI/rigid_flap/OpenFOAM-RigidBody/overview.png differ diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/plotDisplacement.sh b/FSI/rigid_flap/OpenFOAM-RigidBody/plotDisplacement.sh new file mode 100755 index 000000000..894783d70 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/plotDisplacement.sh @@ -0,0 +1,9 @@ +#! /bin/bash +gnuplot -p << EOF +set grid +set title 'Displacement of the rigid body tip' +set xlabel 'time [s]' +set ylabel 'y-displacement [m]' +plot "precice-Solid-watchpoint-flap-tip.log" using 1:5 notitle with lines +EOF + diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/precice-config.xml b/FSI/rigid_flap/OpenFOAM-RigidBody/precice-config.xml new file mode 100644 index 000000000..0d8b5d28f --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/precice-config.xml @@ -0,0 +1,73 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/removeObsoleteFolders.sh b/FSI/rigid_flap/OpenFOAM-RigidBody/removeObsoleteFolders.sh new file mode 100755 index 000000000..6866c8c42 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/removeObsoleteFolders.sh @@ -0,0 +1,22 @@ +#! /bin/bash + +echo "Looking for any time directories without results (e.g. stray functionObjectProperties files, see issue #26 on GitHub)..." + +cd Fluid +for f in [0-9]* [0-9]*.[0-9]*; do + if ! [ -f $f/U ] && ! [ -f $f/T ]; then + rm -rfv $f + fi +done +if [ -d processor0 ]; then + for g in processor*; do + cd $g + for f in [0-9]* [0-9]*.[0-9]*; do + if ! [ -f $f/U ] && ! [ -f $f/T ]; then + rm -rfv $f + fi + done + cd .. + done +fi +cd .. diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/runFluid b/FSI/rigid_flap/OpenFOAM-RigidBody/runFluid new file mode 100755 index 000000000..3223006a0 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/runFluid @@ -0,0 +1,42 @@ +#!/bin/bash +cd ${0%/*} || exit 1 # Run from this directory +. $WM_PROJECT_DIR/bin/tools/RunFunctions # Tutorial run functions + +# Fluid participant + +# Run this script in one terminal and the execute the coupled_elasto_dynamics bindary +# in another terminal. +# These scripts present how the two participants would be started manually. + +# Run this script with "-parallel" for parallel simulations + +# The script "Allclean" cleans-up the result and log files. + +# 1 for true, 0 for false +parallel=0 +if [ "$1" = "-parallel" ]; then + parallel=1 +fi + +echo "Preparing and running the Fluid participant..." + +rm -rfv Fluid/0/ +cp -r Fluid/0.orig/ Fluid/0/ +blockMesh -case Fluid +checkMesh -case Fluid + +# Run +cd Fluid + solver=$(getApplication) + procs=$(getNumberOfProcessors) +cd .. +if [ $parallel -eq 1 ]; then + decomposePar -force -case Fluid + mpirun -np $procs $solver -parallel -case Fluid + reconstructPar -case Fluid +else + $solver -case Fluid +fi + +# Workaround for issue #26 (OF-adapter, relevant for OF .com versions) +./removeObsoleteFolders.sh diff --git a/FSI/rigid_flap/OpenFOAM-RigidBody/runSolid b/FSI/rigid_flap/OpenFOAM-RigidBody/runSolid new file mode 100755 index 000000000..c250d0da5 --- /dev/null +++ b/FSI/rigid_flap/OpenFOAM-RigidBody/runSolid @@ -0,0 +1,10 @@ +#!/bin/bash +cd ${0%/*} || exit 1 # Run from this directory + +# Solid participant + +# Run this script in one terminal and the execute the Fluid participant +# in another terminal. +# These scripts present how the two participants would be started manually. + +./Solid/rigid_body_solver