diff --git a/channel-transport-particles/fluid-nutils/README.md b/channel-transport-particles/fluid-nutils/README.md new file mode 100644 index 000000000..d97cbd5fe --- /dev/null +++ b/channel-transport-particles/fluid-nutils/README.md @@ -0,0 +1,26 @@ +# Hints and TODOs + +## General hints + +- we already have a unidirectional mesh-particle tracing [tutorial](https://precice.org/tutorials-channel-transport-particles.html). This scenario is also available in the [tutorials repository](https://github.com/precice/tutorials/tree/develop/channel-transport-particles) and couples Nutils with MercuryDPM +- the unidirectional tutorial only couples via the fluid velocity and should run out-of-the-box, use `setup-mercurydpm.sh` to run and compile the MercuryDPM participant +- the tutorial uses OpenFOAM or Nutils, and MercuryDPM +- all Nutils simulations use Nutils version 7 +- the bidirectional scenario builds upon this tutorial case + +## Before starting + +- make sure you have the latest develop of preCICE + +## The bidirectional coupling + +- the bidirectional coupling is pull request [706](https://github.com/precice/tutorials/pull/706) +- from the [tutorials](https://github.com/precice/tutorials/tree/develop/channel-bidirectional), checkout the branch `channel-bidirectional` and navigate to `~/tutorials/channel-transport-particles/` +- run the `setup-mercurydpm.sh` script again (it checks out another branch for the solver) +- the overall setup is no too much physics-motivated, using the same boundary conditions leads to the following + +![snapshot-final](../images/tutorials-channel-transport-particles-t-final.png) + +- after you completed the case, I would validate it with a more complex scenario +- while the simulation here is 2D, the final validation case will be 3D +- I put all further comments directly into the `fluid.py` code diff --git a/channel-transport-particles/fluid-nutils/fluid.py b/channel-transport-particles/fluid-nutils/fluid.py index aa91605df..0a27da702 100644 --- a/channel-transport-particles/fluid-nutils/fluid.py +++ b/channel-transport-particles/fluid-nutils/fluid.py @@ -13,9 +13,9 @@ def main(): + # COMMENT: nothing to change here print("Running Nutils") - # define the Nutils mesh nx = 48 ny = 16 step_start = nx // 3 @@ -28,13 +28,19 @@ def main(): step_start:step_end, :step_hight ].withboundary(wall="left,top,right") - # cloud of Gauss points gauss = domain.sample("gauss", degree=4) - # Nutils namespace ns = function.Namespace() ns.x = geom + # COMMENT: looks like the weak is being defined here: the weak form now needs to incorporate the volume fraction + # Note that preCICE/MercuryDPM provides a field called "Alpha", which is the solid volume fraction. The formula + # 7.1 and 7.2 from the PDF use the fluid volume frection (epsilon_f), which is epsilon_f = (1 - Alpha). + # Note that reading from preCICE makes only sense after the initialize call. + # Moreover, we need to incorporate gravity (ε_f*ρ_f*g) in 7.2 and the drag force. The drag force is here zero either + # Note that the drag force is called "ExplicitMomentum" in preCICE. I + # would maybe later still try to implement the explicit-implicit split, + # that's the only reason ns.ubasis = domain.basis("std", degree=2).vector(2) ns.pbasis = domain.basis("std", degree=1) ns.u_i = "ubasis_ni ?u_n" # solution @@ -48,7 +54,7 @@ def main(): ures = gauss.integral("ubasis_ni,j σ_ij d:x" @ ns) pres = gauss.integral("pbasis_n u_k,k d:x" @ ns) - # define Dirichlet boundary condition + # Dirichlet boundary condition sqr = domain.boundary["inflow"].integral("(u_0 - uin)^2 d:x" @ ns, degree=2) sqr += domain.boundary["inflow,outflow"].integral("u_1^2 d:x" @ ns, degree=2) sqr += domain.boundary["wall"].integral("u_k u_k d:x" @ ns, degree=2) @@ -57,13 +63,16 @@ def main(): # preCICE setup participant = precice.Participant("Fluid", "../precice-config.xml", 0, 1) - # define coupling mesh + # define coupling mesh (nothing to change) mesh_name = "Fluid-Mesh" vertices = gauss.eval(ns.x) vertex_ids = participant.set_mesh_vertices(mesh_name, vertices) - # coupling data + # COMMENT: names of the new coupling data fields in the preCICE configuration. The velocity was already there before + solid_fraction_name = "Alpha" + drag_force_name = "ExplicitMomentum" data_name = "Velocity" + pressure_gradient_name = "PressureGradientFull" participant.initialize() @@ -90,13 +99,24 @@ def main(): # potentially adjust non-matching timestep sizes dt = min(solver_dt, precice_dt) + # COMMENT: here we would need to read the solid volume fraction and the drag force from preCICE and incorporate them into the weak form. The pressure gradient would also need to be written to preCICE, but let's start with the velocity first + # Note that the read_data call takes a dt argument, which is the relative read time. For Euler backwards, it would simply be dt + # solid_fraction_values = participant.read_data(mesh_name, solid_fraction_name, vertex_ids, ...) + # drag_force_name = participant.read_data(mesh_name, drag_force_name, vertex_ids, ...) + # Checkpointing for implicit coupling is generally not required + # solve Nutils timestep state["u0"] = state["u"] state["dt"] = dt state = solver.newton(("u", "p"), (ures, pres), constrain=cons, arguments=state).solve(1e-10) velocity_values = gauss.eval(ns.u, **state) + # COMMENT: here we would also need to write the pressure gradient to preCICE + # Important note: OpenFOAM scales the pressue (and thereby also its gradient) by the fluid density, so we would + # need to do the same scaling here. I set the fluid density to 1.2 participant.write_data(mesh_name, data_name, vertex_ids, velocity_values) + # rho_f = 1.2 + # participant.write_data(mesh_name, pressure_gradient_name, vertex_ids, pressure_gradient_values) # do the coupling participant.advance(dt) diff --git a/channel-transport-particles/fluid-openfoam/constant/transportProperties b/channel-transport-particles/fluid-openfoam/constant/transportProperties index 5383adaad..5460c5509 100644 --- a/channel-transport-particles/fluid-openfoam/constant/transportProperties +++ b/channel-transport-particles/fluid-openfoam/constant/transportProperties @@ -8,4 +8,7 @@ FoamFile transportModel Newtonian; -nu nu [ 0 2 -1 0 0 0 0 ] 1; +nu [0 2 -1 0 0 0 0] 1.5e-05; // = air at 20°C m^2/s +rho 1.2; // kg/m^3 +g (0 -9.81 0); + diff --git a/channel-transport-particles/fluid-openfoam/system/controlDict b/channel-transport-particles/fluid-openfoam/system/controlDict index 55c798b99..88734cfbf 100644 --- a/channel-transport-particles/fluid-openfoam/system/controlDict +++ b/channel-transport-particles/fluid-openfoam/system/controlDict @@ -6,7 +6,7 @@ FoamFile object controlDict; } -application pimpleFoam; +application AndersonJacksonFoam; startFrom startTime; diff --git a/channel-transport-particles/fluid-openfoam/system/fvSchemes b/channel-transport-particles/fluid-openfoam/system/fvSchemes index 80c096192..fc2e29bcd 100644 --- a/channel-transport-particles/fluid-openfoam/system/fvSchemes +++ b/channel-transport-particles/fluid-openfoam/system/fvSchemes @@ -18,7 +18,7 @@ gradSchemes divSchemes { - default none; + default Gauss linear; div(phi,U) bounded Gauss upwind; div((nuEff*dev2(T(grad(U))))) Gauss linear; } diff --git a/channel-transport-particles/fluid-openfoam/system/preciceDict b/channel-transport-particles/fluid-openfoam/system/preciceDict index 7368a7a56..506e26435 100644 --- a/channel-transport-particles/fluid-openfoam/system/preciceDict +++ b/channel-transport-particles/fluid-openfoam/system/preciceDict @@ -12,6 +12,14 @@ participant Fluid; modules (FF); +FF +{ + nameAlpha alphaSolid; + nameExplicitMomentum ExplicitMomentum; + nameImplicitMomentum ImplicitMomentum; +}; + + interfaces { Interface1 @@ -22,11 +30,15 @@ interfaces readData ( + ImplicitMomentum + ExplicitMomentum + Alpha ); writeData ( Velocity + PressureGradientFull ); }; -}; \ No newline at end of file +}; diff --git a/channel-transport-particles/images/tutorials-channel-transport-particles-t-final.png b/channel-transport-particles/images/tutorials-channel-transport-particles-t-final.png new file mode 100644 index 000000000..a6abbdad2 Binary files /dev/null and b/channel-transport-particles/images/tutorials-channel-transport-particles-t-final.png differ diff --git a/channel-transport-particles/particles-mercurydpm/setup-mercurydpm.sh b/channel-transport-particles/particles-mercurydpm/setup-mercurydpm.sh index 1a5863861..4d04d4a25 100644 --- a/channel-transport-particles/particles-mercurydpm/setup-mercurydpm.sh +++ b/channel-transport-particles/particles-mercurydpm/setup-mercurydpm.sh @@ -4,7 +4,7 @@ set -euo pipefail # ---- Config (you can override via env vars) ---- REPO_HTTPS_DEFAULT="https://bitbucket.org/davidscn/mercurydpm.git" REPO_URL="${REPO_URL:-$REPO_HTTPS_DEFAULT}" # or git@bitbucket.org:davidscn/mercurydpm.git -BRANCH="${BRANCH:-channel-transport-tutorial}" +BRANCH="${BRANCH:-channel-bidirectional}" SRC_DIR="${SRC_DIR:-$PWD/mercurydpm}" BUILD_DIR="${BUILD_DIR:-$SRC_DIR/build}" diff --git a/channel-transport-particles/precice-config.xml b/channel-transport-particles/precice-config.xml index 01af32b7a..aff17a43c 100644 --- a/channel-transport-particles/precice-config.xml +++ b/channel-transport-particles/precice-config.xml @@ -1,43 +1,84 @@ - + + filter="(%Severity% > debug) and not ((%Rank% != 0))" /> - + - + + + + + + + + + + + + + - + + + + + - + + - + + - - - + + + + + +