From c427c1f7d362b6c8ce4712d25237a5cad219b676 Mon Sep 17 00:00:00 2001
From: max <maxwesterlund2@gmail.com>
Date: Sun, 3 May 2026 02:10:28 +0200
Subject: [PATCH] Added pipe friction

---
 Components/Pipe1D.cs | 26 ++++++++++++++++++++------
 Core/Simulation.cs   |  8 +++++---
 2 files changed, 25 insertions(+), 9 deletions(-)

diff --git a/Components/Pipe1D.cs b/Components/Pipe1D.cs
index fd19cb3..dea0d24 100644
--- a/Components/Pipe1D.cs
+++ b/Components/Pipe1D.cs
@@ -15,9 +15,10 @@ namespace FluidSim.Components
         public Port PortA { get; }
         public Port PortB { get; }
         public double Area => _area;
+        public double DampingMultiplier { get; set; } = 1.0;
 
         private int _n;
-        private double _dx, _dt, _gamma, _area;
+        private double _dx, _dt, _gamma, _area, _diameter;
         private double[] _rho, _rhou, _E;
 
         private double _rhoLeft, _pLeft;
@@ -33,8 +34,6 @@ namespace FluidSim.Components
         private const double CflTarget = 0.8;
         private const double ReferenceSoundSpeed = 340.0;
 
-        public double FrictionFactor { get; set; } = 0.02;
-
         public int GetCellCount() => _n;
         public double GetCellDensity(int i) => _rho[i];
         public double GetCellPressure(int i) => Pressure(i);
@@ -66,6 +65,9 @@ namespace FluidSim.Components
             _area = area;
             _gamma = 1.4;
 
+            // Hydraulic diameter for a circular pipe
+            _diameter = 2.0 * Math.Sqrt(area / Math.PI);
+
             _rho = new double[_n];
             _rhou = new double[_n];
             _E = new double[_n];
@@ -229,9 +231,14 @@ namespace FluidSim.Components
                     break;
             }
 
-            // Cell update
+            // ---- Cell update with linear laminar damping ----
+            double radius = _diameter / 2.0;
+            double mu_air = 1.8e-5;                       // dynamic viscosity of air (Pa·s)
+            double laminarCoeff = DampingMultiplier * 8.0 * mu_air / (radius * radius);
+
             for (int i = 0; i < n; i++)
             {
+                // Flux divergence
                 double dM = (Fm[i + 1] - Fm[i]) / _dx;
                 double dP = (Fp[i + 1] - Fp[i]) / _dx;
                 double dE = (Fe[i + 1] - Fe[i]) / _dx;
@@ -240,15 +247,22 @@ namespace FluidSim.Components
                 _rhou[i] -= dtSub * dP;
                 _E[i] -= dtSub * dE;
 
-                if (_rho[i] < 1e-12) _rho[i] = 1e-12;
+                // Laminar viscous damping on momentum (implicit exponential decay)
+                double rho = Math.Max(_rho[i], 1e-12);
+                double dampingRate = laminarCoeff / rho;   // 1/s
+                double dampingFactor = Math.Exp(-dampingRate * dtSub);
+                _rhou[i] *= dampingFactor;
+                // Note: total energy _E[i] is unchanged – kinetic energy loss becomes internal heat
 
+                // Physical bounds
+                if (_rho[i] < 1e-12) _rho[i] = 1e-12;
                 double kinetic = 0.5 * _rhou[i] * _rhou[i] / _rho[i];
                 double pMin = 100.0;
                 double eMin = pMin / ((_gamma - 1) * _rho[i]) + kinetic;
                 if (_E[i] < eMin) _E[i] = eMin;
             }
 
-            // Port quantities (only meaningful for volume coupled ends)
+            // Port quantities (only meaningful for volume‑coupled ends)
             double mdotA_sub = _leftBCType == BoundaryType.VolumeCoupling && _leftBCSet ? Fm[0] * _area : 0.0;
             double mdotB_sub = _rightBCType == BoundaryType.VolumeCoupling && _rightBCSet ? -Fm[n] * _area : 0.0;
 
diff --git a/Core/Simulation.cs b/Core/Simulation.cs
index 308f058..9188aed 100644
--- a/Core/Simulation.cs
+++ b/Core/Simulation.cs
@@ -15,17 +15,18 @@ namespace FluidSim.Core
         private static float sample;
         private static double ambientPressure = 1.0 * Units.atm;
 
+        private static bool enableLogging = false;
+
         public static void Initialize(int sampleRate)
         {
             dt = 1.0 / sampleRate;
 
-            double length = 0.2;
-            double radius = 5 * Units.mm;
+            double length = 2;
+            double radius = 20 * Units.mm;
             double area = Units.AreaFromDiameter(radius);
 
             pipe = new Pipe1D(length, area, sampleRate, forcedCellCount: 80);
             pipe.SetUniformState(1.225, 0.0, ambientPressure);
-            pipe.FrictionFactor = 0.0;
 
             solver = new Solver();
             solver.SetTimeStep(dt);
@@ -56,6 +57,7 @@ namespace FluidSim.Core
 
         public static void Log()
         {
+            if (!enableLogging) return;
             if (stepCount % 10 == 0 && stepCount < 1000)
             {
                 double pMid = pipe.GetPressureAtFraction(0.5);