Water Hammer Tutorial

water-hammer/README.md

@@ -0,0 +1,166 @@
+---
+title: Partitioned Water Hammer
+keywords: OpenFOAM, Nutils, preCICE, multiscale, fluid, transient
+summary: The Partitioned Water Hammer tutorial simulates unsteady pressure wave propagation in pipe systems using different 1D and 3D configurations coupled via preCICE.
+---
+
+## Setup
+
+We solve a **partitioned water hammer problem** using a 1D–3D coupling approach.  
+In this tutorial, the computational domain is split into two coupled regions: a 1D pipe section and a 3D pipe section.  
+The coupling is performed using **preCICE**.
+
+Case setup inspired by [1].  
+This tutorial extends the study conducted in [2], which implemented the water hammer benchmark using a 1D-3D coupling with preCICE.  
+In that study, the cross-section of the pipe was squared. It has been changed to a circular cross-section in the present tutorial.  
+`1D` denotes the reduced-order domain (e.g., a Nutils solver) and `3D` denotes the full 3D CFD domain (e.g., OpenFOAM).
+
+The problem consists of a straight pipe of length `L = 1000 m` and diameter `D = 2 m`.  
+We partition the domain at `y_c = 500 m`, where the coupling interface is located.  
+The **1D domain** solves the unsteady compressible flow equations using Nutils, while the **3D domain** is solved using OpenFOAM.  
+Both solvers are coupled via preCICE by exchanging the **pressure** and **axial velocity** at the interface.
+
+Two coupling directions are possible:  
+- **1D → 3D**: The 1D solver provides the interface pressure to the 3D solver, which responds with velocity.  
+- **3D → 1D**: The 3D solver provides the pressure, and the 1D solver returns the velocity.
+
+The initial inlet pressure is set to `p_in = 98100 Pa`.  
+The opening valve at the outlet is modeled through a prescribed outlet velocity, which increases linearly from `0` to `1 m/s` over the first `5 s` and remains constant afterwards.  
+This sudden valve opening generates pressure disturbances that propagate through the pipe, resulting in the characteristic **pressure wave oscillations** known as the *water hammer* phenomenon.
+
+## Configuration
+
+preCICE configuration for the 1D-3D simulation (image generated using the [precice-config-visualizer](https://precice.org/tooling-config-visualization.html))
+
+![preCICE configuration visualization 1D-3D](images/tutorials-waterhammer-1d3d-precice-config.png)
+
+preCICE configuration for the 3D-1D simulation (image generated using the [precice-config-visualizer](https://precice.org/tooling-config-visualization.html))
+
+![preCICE configuration visualization 3D-1D](images/tutorials-waterhammer-3d1d-precice-config.png)
+
+## Available solvers
+
+* OpenFOAM (sonicLiquidFoam). A compressible OpenFOAM solver. For more information, have a look at the [OpenFOAM adapter documentation](https://precice.org/adapter-openfoam-overview.html)
+* Nutils v9. A Python-based finite element framework. For more information, see the [Nutils adapter documentation](https://precice.org/adapter-nutils.html)
+
+## Running the Simulation
+
+First, select which coupling you want to run. This sets the correct `precice-config.xml` symlink:
+
+```bash
+# Choose one (1D → 3D or 3D → 1D)
+./setcase.sh 1d3d
+# or
+./setcase.sh 3d1d
+```
+
+Open **two terminals** and start the corresponding participants for your chosen setup.
+
+### Example A — 1D → 3D coupling
+
+Terminal 1:
+```bash
+cd fluid1d-left
+./run.sh
+```
+
+Terminal 2:
+```bash
+cd fluid3d-right
+./run.sh
+```
+
+### Example B — 3D → 1D coupling
+
+Terminal 1:
+```bash
+cd fluid3d-left
+./run.sh
+```
+
+Terminal 2:
+```bash
+cd fluid1d-right
+./run.sh
+```
+
+> Tip: If you switch coupling direction later, rerun `./setcase.sh` with the other option before launching the participants.
+
+## Visualization
+
+The output of the coupled simulation is written into the folders `fluid1d-left`, `fluid1d-right`, `fluid3d-left`, and `fluid3d-right`, depending on which coupling direction (`1d3d` or `3d1d`) you selected.
+
+### 3D domain (OpenFOAM)
+
+For the 3D participant, all simulation results are stored in the time directories inside the respective case folder (e.g., `fluid3d-right/`).  
+You can visualize the flow field and pressure distribution using **ParaView** by opening the case file:
+
+```bash
+paraview fluid3d-right/fluid3d-right.foam
+```
+
+or, for the left domain if applicable:
+
+```bash
+paraview fluid3d-left/fluid3d-left.foam
+```
+
+Typical fields to inspect include:  
+- `p` – pressure (to observe the pressure wave propagation)  
+- `U` – velocity (to visualize the flow development near the coupling interface)
+
+In ParaView, you can also create line plots along the pipe centerline or use temporal filters to observe wave propagation over time.
+
+We also record pressure and velocity at fixed points each time step using the OpenFOAM `probes` function object.
+
+**Probe setup (excerpt):**
+```c
+#includeEtc "caseDicts/postProcessing/probes/probes.cfg"
+
+fields (p U);
+probeLocations
+(
+    (0 999 0)
+    (0 500 0)
+);
+// For the left 3D domain use instead:
+// probeLocations ((0 0 0) (0 500 0));
+```
+
+**Output location:**  
+- `fluid3d-right/postProcessing/probes/0/p`  
+- `fluid3d-right/postProcessing/probes/0/U`
+
+Each file contains a header (commented with `#`) and time-series columns for each probe.
+
+### 1D domain (Nutils)
+
+The 1D participant writes results to a file named `watchpoint.txt`, containing the temporal evolution of key quantities at selected spatial locations.
+
+The structure of each line in the file is:
+
+```
+time; p_in; u_in; p_out; u_out; p_mid; u_mid
+```
+
+where:  
+- `p_in`, `u_in` → pressure and velocity at the inlet of the 1D domain,  
+- `p_out`, `u_out` → pressure and velocity at the outlet of the 1D domain,  
+- `p_mid`, `u_mid` → pressure and velocity at the midpoint of the 1D domain.
+
+### Example visualization
+
+![Pressure evolution at the outlet of the 3D domain in the 1D-3D simulation](images/tutorials-waterhammer-1d3d-outlet-pressure.png)
+
+Pressure evolution at the outlet of the 3D domain during the 1D–3D water hammer simulation.  
+
+## References
+
+[1] C. Wang, H. Nilsson, J. Yang, *et al.*  
+**1D–3D coupling for hydraulic system transient simulations.**  
+*Computer Physics Communications*, 210:1–9, 2017. 
+[https://doi.org/10.1016/j.cpc.2016.09.007](https://doi.org/10.1016/j.cpc.2016.09.007) 
+
+[2] G. Chourdakis, B. Uekermann, G. van Zwieten, and H. van Brummelen.  
+**Coupling OpenFOAM to different solvers, physics, models, and dimensions using preCICE.**  
+[https://mediatum.ub.tum.de/doc/1515271/1515271.pdf](https://mediatum.ub.tum.de/doc/1515271/1515271.pdf)

water-hammer/clean-tutorial.sh

@@ -0,0 +1,11 @@
+#!/usr/bin/env sh
+set -e -u
+
+# shellcheck disable=SC1091
+. ../tools/cleaning-tools.sh
+
+clean_tutorial .
+clean_precice_logs .
+rm -fv ./*.log
+rm -fv ./*.vtu
+

water-hammer/fluid1d-left/Fluid1D.py

@@ -0,0 +1,147 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import treelog
+from nutils import mesh, function, solver, cli
+import precice
+
+
+def main(nelems=200, dt=.005, refdensity=1e3, refpressure=101325.0, psi=1e-6, viscosity=1e-3, theta=0.5):
+
+    # --- preCICE initialization ---
+
+    participant_name = "Fluid1DLeft"
+    config_file_name = "../precice-config.xml"
+    solver_process_index = 0
+    solver_process_size = 1
+    participant = precice.Participant(participant_name, config_file_name, solver_process_index, solver_process_size)
+    mesh_name = "Fluid1DLeft-Mesh"
+    velocity_name = "Velocity"
+    pressure_name = "Pressure"
+    positions = [[0.0, 500.0, 0.0]]     # Define a single coupling point (3D format required by preCICE)
+    vertex_ids = participant.set_mesh_vertices(mesh_name, positions)
+
+    participant.initialize()
+    precice_dt = participant.get_max_time_step_size()
+
+    # --- Nutils domain setup ---
+
+    domain, geom = mesh.rectilinear([np.linspace(0, 500, nelems + 1)])    # 1D domain from 0 to 500 with nelems elements
+
+    ns = function.Namespace()
+    ns.x = geom
+    ns.dt = dt
+    ns.θ = theta
+    ns.ρref = refdensity
+    ns.pref = refpressure
+    ns.pin = 98100  # Inlet pressure (Pa)
+    ns.μ = viscosity
+    ns.ψ = psi
+
+    # Define basis functions: velocity (vector, degree 2), density (scalar, degree 1)
+    ns.ubasis, ns.ρbasis = function.chain([
+        domain.basis('std', degree=2).vector(domain.ndims),
+        domain.basis('std', degree=1)
+    ])
+
+    # Define trial and previous-step functions
+    ns.u_i = 'ubasis_ni ?lhs_n'
+    ns.u0_i = 'ubasis_ni ?lhs0_n'
+    ns.ρ = 'ρref + ρbasis_n ?lhs_n'
+    ns.ρ0 = 'ρref + ρbasis_n ?lhs0_n'
+    ns.p = 'pref + (ρ - ρref) / ψ'
+    ns.utheta_i = 'θ u_i + (1 - θ) u0_i'
+    ns.ρtheta = 'θ ρ + (1 - θ) ρ0'
+
+    # Cauchy stress tensor: viscous + pressure
+    ns.σ_ij = 'μ (utheta_i,j + utheta_j,i) - p δ_ij'
+
+    # --- Weak form residuals ---
+
+    # Mass conservation
+    res = domain.integral(
+        'ρbasis_n ((ρ - ρ0) / dt + (ρtheta utheta_k)_,k) d:x' @ ns,
+        degree=4
+    )
+
+    # Momentum conservation: unsteady + convective + diffusive terms
+    res += domain.integral(
+        '(ubasis_ni (((ρ u_i - ρ0 u0_i) / dt) + (ρtheta utheta_i utheta_j)_,j) + ubasis_ni,j σ_ij) d:x' @ ns,
+        degree=4
+    )
+
+    # Weakly enforced inlet pressure boundary condition
+    res += domain.boundary['left'].integral(
+        'pin ubasis_ni n_i d:x' @ ns,
+        degree=4
+    )
+
+    # Initial condition
+    lhs0 = np.zeros(res.shape)
+
+    # Time-stepping
+    t = 0.0
+    timestep = 0
+    bezier = domain.sample('bezier', 2)
+
+    f = open("watchpoint.txt", "w")
+
+    # --- Time loop with preCICE coupling ---
+    while participant.is_coupling_ongoing():
+
+        # Read velocity data from other solver via preCICE
+        u_read = participant.read_data(mesh_name, velocity_name, vertex_ids, precice_dt)
+
+        # Constrain outlet velocity to match coupled value
+        stringintegral = f'(u_0 - ({u_read[0][1]}))^2 d:x'
+        sqr = domain.boundary['right'].integral(stringintegral @ ns, degree=4)
+        cons = solver.optimize('lhs', sqr, droptol=1e-14)
+
+        # Write checkpoint if required by preCICE
+        if participant.requires_writing_checkpoint():
+            lhscheckpoint = lhs0
+
+        # Solve nonlinear system (Newton's method)
+        with treelog.context('stokes'):
+            lhs = solver.newton(
+                'lhs',
+                residual=res,
+                constrain=cons,
+                arguments=dict(lhs0=lhs0),
+                lhs0=lhs0
+            ).solve(1e-08)
+
+            # Evaluate fields for visualization/debugging
+            x, p, u, ρ = bezier.eval(['x_i', 'p', 'u_i', 'ρ'] @ ns, arguments=dict(lhs=lhs))
+
+        # Send pressure at the right boundary to the other solver
+        write_press = [[p[-1]]]
+        participant.write_data(mesh_name, pressure_name, vertex_ids, write_press)
+
+        # Advance in pseudo-time
+        participant.advance(precice_dt)
+        precice_dt = participant.get_max_time_step_size()
+
+        # Restore checkpoint if needed
+        if participant.requires_reading_checkpoint():
+            lhs0 = lhscheckpoint
+        else:
+            # Update old solution
+            lhs0 = lhs
+            timestep += timestep
+
+            # Save probe values (time, inlet pressure, inlet velocity, outlet
+            # pressure, outlet velocity, pressure at the middle, velocity at the
+            # middle)
+            x, p, ρ, u = bezier.eval(['x_i', 'p', 'ρ', 'u_i'] @ ns, lhs=lhs)
+            f.write("%e; %e; %e; %e; %e; %e; %e\n" % (t, p[0], u[0], p[-1], u[-1], p[199], u[199]))
+            f.flush()
+
+            t += precice_dt  # advance pseudo-time
+
+    # Finalize preCICE
+    participant.finalize()
+    f.close()
+
+
+if __name__ == '__main__':
+    cli.run(main)

water-hammer/fluid1d-left/clean.sh

@@ -0,0 +1,8 @@
+#!/usr/bin/env sh
+set -e -u
+
+. ../../tools/cleaning-tools.sh
+
+rm -f ./results/Fluid1D_*
+clean_precice_logs .
+clean_case_logs .

water-hammer/fluid1d-left/requirements.txt

@@ -0,0 +1,6 @@
+setuptools # required by nutils
+nutils==9
+numpy >1, <2
+pyprecice~=3.0
+matplotlib
+treelog

water-hammer/fluid1d-left/run.sh

@@ -0,0 +1,13 @@
+#!/usr/bin/env bash
+set -e -u
+
+. ../../tools/log.sh
+exec > >(tee --append "$LOGFILE") 2>&1
+
+python3 -m venv .venv
+. .venv/bin/activate
+pip install -r requirements.txt && pip freeze > pip-installed-packages.log
+
+NUTILS_RICHOUTPUT=no python3 Fluid1D.py
+
+close_log