diff --git a/CtrEditor.csproj b/CtrEditor.csproj index 5593f64..b8c6be5 100644 --- a/CtrEditor.csproj +++ b/CtrEditor.csproj @@ -188,6 +188,7 @@ + diff --git a/Documentation/HydraulicSimulator DOC.md b/Documentation/HydraulicSimulator DOC.md new file mode 100644 index 0000000..e64de8a --- /dev/null +++ b/Documentation/HydraulicSimulator DOC.md @@ -0,0 +1,397 @@ +# 📖 Documentación - Simulador Hidráulico C# .NET 8 + +## 🎯 Descripción General + +Esta librería permite simular redes hidráulicas complejas con bombas, tuberías, válvulas y tanques. Está diseñada para calcular flujos y presiones en cada punto de la red, permitiendo posteriormente calcular niveles en tanques basados en los flujos calculados. + +## 🏗️ Arquitectura de la Librería + +### 📁 Estructura del Proyecto +``` +HydraulicSimulator/ +├── Models/ # Modelos principales +│ ├── Fluid.cs # Propiedades del fluido +│ ├── Node.cs # Nodos de la red +│ ├── Branch.cs # Ramas que conectan nodos +│ ├── Element.cs # Clase base para elementos +│ ├── Pipe.cs # Tuberías +│ ├── PumpHQ.cs # Bombas con curva H-Q +│ ├── ValveKv.cs # Válvulas con coeficiente Kv +│ ├── MinorLoss.cs # Pérdidas menores +│ └── HydraulicNetwork.cs # Red hidráulica y solver +├── Core/ # Funcionalidades auxiliares +│ ├── TestRunner.cs # Casos de prueba +│ └── InteractiveSimulator.cs # Simulador interactivo +└── Program.cs # Punto de entrada +``` + +## 🧱 Componentes Principales + +### 1. **Fluid** - Propiedades del Fluido +```csharp +var fluid = new Fluid( + rho: 1000, // Densidad (kg/m³) + mu: 0.001, // Viscosidad dinámica (Pa·s) + name: "Water", // Nombre del fluido + temp: 20.0, // Temperatura (°C) + pVapor: 2337.0 // Presión de vapor (Pa) +); + +// Fluido predefinido +var water = Fluid.Water20C; // Agua a 20°C +``` + +### 2. **Node** - Nodos de la Red +```csharp +// Nodo con presión fija (tanques, descargas) +var tank = new Node("TANQUE_01", fixedP: true, p: 0); + +// Nodo con presión libre (calculada por el solver) +var junction = new Node("UNION_01", fixedP: false, p: 0); +``` + +**Propiedades:** +- `Name`: Nombre identificador +- `FixedP`: Si la presión es fija (true) o calculada (false) +- `P`: Presión en Pascal (Pa) + +### 3. **Pipe** - Tuberías +```csharp +var pipe = new Pipe( + length: 100, // Longitud (m) + diameter: 0.1, // Diámetro interno (m) + roughness: 0.0015 // Rugosidad absoluta (m) +); +``` + +**Características:** +- Usa ecuación de Darcy-Weisbach +- Factor de fricción por Swamee-Jain +- Cálculo automático de pérdidas por fricción + +### 4. **PumpHQ** - Bombas +```csharp +var pump = new PumpHQ( + h0: 100, // Cabeza a caudal cero (m) + q0: 0.02, // Caudal a cabeza cero (m³/s) + speedRel: 1.0, // Velocidad relativa (1.0 = nominal) + direction: 1 // Dirección (+1 o -1) +); +``` + +**Curva característica:** `H = H0 * (1 - (Q/Q0)²)` + +### 5. **ValveKv** - Válvulas +```csharp +var valve = new ValveKv( + kvFull: 25, // Kv a apertura completa (m³/h/√bar) + opening: 1.0 // Apertura (0.0 a 1.0) +); +``` + +### 6. **MinorLoss** - Pérdidas Menores +```csharp +var elbow = new MinorLoss(k: 0.9); // Codo 90° +var tee = new MinorLoss(k: 1.8); // Té +``` + +## 🔧 Uso Básico en tu Aplicación + +### 1. **Crear una Red Hidráulica** +```csharp +using HydraulicSimulator.Models; + +// 1. Crear la red +var network = new HydraulicNetwork(); + +// 2. Agregar nodos +network.AddNode("TANQUE_A", 0); // Tanque a presión atmosférica +network.AddNode("UNION_1"); // Nodo intermedio +network.AddNode("PROCESO_1", 50000); // Proceso a 0.5 bar +network.AddNode("DESCARGA", 0); // Descarga a atmósfera + +// 3. Crear elementos +var pump = new PumpHQ(h0: 80, q0: 0.01); +var mainPipe = new Pipe(length: 50, diameter: 0.08, roughness: 0.0015); +var valve = new ValveKv(kvFull: 20, opening: 0.8); +var dischargePipe = new Pipe(length: 30, diameter: 0.06, roughness: 0.0015); + +// 4. Conectar elementos en ramas +network.AddBranch("TANQUE_A", "UNION_1", + new List { pump, mainPipe }); + +network.AddBranch("UNION_1", "PROCESO_1", + new List { valve }); + +network.AddBranch("PROCESO_1", "DESCARGA", + new List { dischargePipe }); + +// 5. Resolver la red +var result = network.Solve(); + +if (result.Converged) +{ + Console.WriteLine("✅ Simulación exitosa"); + network.Report(); // Mostrar resultados +} +``` + +### 2. **Obtener Resultados** +```csharp +// Flujos en cada rama (m³/s) +foreach (var branch in network.Branches) +{ + Console.WriteLine($"Rama {branch.Name}: {branch.Q:F6} m³/s"); +} + +// Presiones en cada nodo (Pa) +foreach (var node in network.Nodes) +{ + double pressureBar = node.Value.P / 100000.0; + Console.WriteLine($"Nodo {node.Key}: {pressureBar:F2} bar"); +} + +// Resultados específicos +var flowRate = result.Flows["TANQUE_A->UNION_1"]; +var pressure = result.Pressures["UNION_1"]; +``` + +## 🎮 Integración con Objetos Gráficos + +### Mapeo de Objetos Gráficos a Elementos +```csharp +public class GraphicToHydraulicMapper +{ + public HydraulicNetwork CreateNetworkFromGraphics(List graphics) + { + var network = new HydraulicNetwork(); + + // 1. Agregar todos los nodos + foreach (var graphic in graphics.Where(g => g.Type == "Node")) + { + var isFixed = graphic.Properties.ContainsKey("FixedPressure"); + var pressure = isFixed ? Convert.ToDouble(graphic.Properties["Pressure"]) : 0; + + network.AddNode(graphic.Id, isFixed ? pressure : null); + } + + // 2. Crear elementos hidráulicos desde gráficos + foreach (var graphic in graphics.Where(g => g.Type != "Node")) + { + Element element = graphic.Type switch + { + "Pump" => new PumpHQ( + h0: GetProperty(graphic, "Head"), + q0: GetProperty(graphic, "MaxFlow") / 3600.0 // Convertir m³/h a m³/s + ), + "Pipe" => new Pipe( + length: GetProperty(graphic, "Length"), + diameter: GetProperty(graphic, "Diameter"), + roughness: GetProperty(graphic, "Roughness", 0.0015) + ), + "Valve" => new ValveKv( + kvFull: GetProperty(graphic, "Kv"), + opening: GetProperty(graphic, "Opening", 1.0) + ), + _ => throw new NotSupportedException($"Elemento {graphic.Type} no soportado") + }; + + // 3. Conectar elemento entre nodos + network.AddBranch( + graphic.FromNodeId, + graphic.ToNodeId, + new List { element }, + graphic.Id + ); + } + + return network; + } + + private T GetProperty(GraphicElement graphic, string key, T defaultValue = default) + { + return graphic.Properties.ContainsKey(key) + ? (T)Convert.ChangeType(graphic.Properties[key], typeof(T)) + : defaultValue; + } +} +``` + +## 🏭 Casos Especiales: Tanques y Descargas + +### 1. **Tanque de Suministro** +```csharp +// Tanque con nivel fijo (presión constante) +network.AddNode("TANQUE_SUMINISTRO", 0); // Presión atmosférica + +// Tanque con altura (presión hidrostática) +double height = 10; // metros +double pressure = 1000 * 9.81 * height; // ρgh +network.AddNode("TANQUE_ELEVADO", pressure); +``` + +### 2. **Tanque de Proceso (Volumen Variable)** +```csharp +// Nodo con presión libre - el nivel se calcula después +network.AddNode("TANQUE_PROCESO"); // Presión calculada por el solver + +// Después de la simulación, calcular nivel +public double CalculateLevel(double pressure, double density = 1000) +{ + return pressure / (density * 9.81); // h = P/(ρg) +} +``` + +### 3. **Descarga a Atmósfera** +```csharp +// Descarga libre (presión atmosférica) +network.AddNode("DESCARGA", 0); + +// Descarga con contrapresión +network.AddNode("DESCARGA_PRESURIZADA", 20000); // 0.2 bar +``` + +### 4. **Tanque con Cálculo de Nivel Dinámico** +```csharp +public class Tank +{ + public string NodeId { get; set; } + public double Area { get; set; } // m² - área del tanque + public double CurrentVolume { get; set; } // m³ - volumen actual + public double Height => CurrentVolume / Area; // m - altura actual + + // Actualizar volumen basado en flujos + public void UpdateVolume(Dictionary flows, double deltaTime) + { + double netFlow = 0; + + // Sumar flujos entrantes, restar salientes + foreach (var flow in flows) + { + if (flow.Key.EndsWith($"->{NodeId}")) + netFlow += flow.Value; // Entrante + else if (flow.Key.StartsWith($"{NodeId}->")) + netFlow -= flow.Value; // Saliente + } + + // Actualizar volumen + CurrentVolume += netFlow * deltaTime; + CurrentVolume = Math.Max(0, CurrentVolume); // No puede ser negativo + } + + // Calcular presión en el fondo del tanque + public double BottomPressure(double density = 1000) + { + return density * 9.81 * Height; // Presión hidrostática + } +} +``` + +## 📊 Ejemplo Completo: Sistema de Mixers + +```csharp +public class MixerSystemSimulation +{ + public void RunSimulation() + { + var network = new HydraulicNetwork(); + + // Nodos del sistema + network.AddNode("TANQUE_PRINCIPAL", 0); // Suministro + network.AddNode("MIXER_1", 100000); // Proceso 1 (1 bar) + network.AddNode("MIXER_2", 100000); // Proceso 2 (1 bar) + network.AddNode("DESCARGA", 50000); // Descarga (0.5 bar) + + // Equipos + var pump1 = new PumpHQ(h0: 120, q0: 0.009); + var pump2 = new PumpHQ(h0: 110, q0: 0.008); + var supply1 = new Pipe(length: 50, diameter: 0.1, roughness: 0.0015); + var supply2 = new Pipe(length: 45, diameter: 0.08, roughness: 0.0015); + var valve1 = new ValveKv(kvFull: 25, opening: 1.0); + var valve2 = new ValveKv(kvFull: 20, opening: 1.0); + var discharge1 = new Pipe(length: 20, diameter: 0.06, roughness: 0.0015); + var discharge2 = new Pipe(length: 25, diameter: 0.06, roughness: 0.0015); + + // Conexiones + network.AddBranch("TANQUE_PRINCIPAL", "MIXER_1", + new List { pump1, supply1, valve1 }); + network.AddBranch("TANQUE_PRINCIPAL", "MIXER_2", + new List { pump2, supply2, valve2 }); + network.AddBranch("MIXER_1", "DESCARGA", + new List { discharge1 }); + network.AddBranch("MIXER_2", "DESCARGA", + new List { discharge2 }); + + // Simular + var result = network.Solve(); + + if (result.Converged) + { + // Procesar resultados + ProcessResults(result); + } + else + { + Console.WriteLine($"❌ No convergió después de {result.Iterations} iteraciones"); + } + } + + private void ProcessResults(SolutionResult result) + { + Console.WriteLine("🎯 RESULTADOS DE SIMULACIÓN:"); + + foreach (var flow in result.Flows) + { + double flowM3h = flow.Value * 3600; // Convertir a m³/h + Console.WriteLine($" {flow.Key}: {flowM3h:F2} m³/h"); + } + + foreach (var pressure in result.Pressures) + { + double pressureBar = pressure.Value / 100000.0; // Convertir a bar + Console.WriteLine($" {pressure.Key}: {pressureBar:F2} bar"); + } + } +} +``` + +## 🚀 Ejecución desde Línea de Comandos + +```bash +# Casos de prueba disponibles +dotnet run test simple-pipe --verbose +dotnet run test mixer-system --json +dotnet run test complex-network --flow 25 --pressure 30 + +# Ayuda +dotnet run help +``` + +## ⚙️ Parámetros de Configuración + +### Solver Parameters +```csharp +var result = network.Solve( + maxIterations: 200, // Máximo de iteraciones + tolerance: 1e-4, // Tolerancia de convergencia + relaxationFactor: 0.7, // Factor de relajación + verbose: true // Mostrar progreso +); +``` + +### Unidades del Sistema +- **Presión**: Pascal (Pa) - Convertir: `bar = Pa / 100000` +- **Flujo**: m³/s - Convertir: `m³/h = m³/s * 3600` +- **Longitud**: metros (m) +- **Diámetro**: metros (m) +- **Rugosidad**: metros (m) + +## 🎯 Recomendaciones para tu Aplicación + +1. **📋 Mapea tus objetos gráficos** a los elementos hidráulicos +2. **🔄 Implementa actualización dinámica** de niveles de tanques +3. **⚡ Usa tolerancias apropiadas** según la complejidad de tu red +4. **📊 Implementa visualización** de resultados en tiempo real +5. **💾 Guarda configuraciones** para reutilizar casos de simulación + +¿Te gustaría que desarrolle alguna sección específica o que agregue ejemplos para casos particulares de tu aplicación? diff --git a/HydraulicSimulator/Models/Branch.cs b/HydraulicSimulator/Models/Branch.cs new file mode 100644 index 0000000..1a6437e --- /dev/null +++ b/HydraulicSimulator/Models/Branch.cs @@ -0,0 +1,82 @@ +using System; +using System.Collections.Generic; +using System.Linq; + +namespace HydraulicSimulator.Models +{ + /// + /// Rama (elementos en serie) + /// + public class Branch + { + public string N1 { get; set; } + public string N2 { get; set; } + public List Elements { get; set; } + public string Name { get; set; } + public double Q { get; set; } = 0.0; // caudal actual (signo n1->n2) + + public Branch(string n1, string n2, List elements, string name = "") + { + N1 = n1; + N2 = n2; + Elements = new List(elements); + Name = string.IsNullOrEmpty(name) ? $"{n1}->{n2}" : name; + } + + public double Dp(double q, Fluid fluid) + { + return Elements.Sum(e => e.Dp(q, fluid)); + } + + public double DdpDq(double q, Fluid fluid) + { + return Elements.Sum(e => e.DdpDq(q, fluid)); + } + + /// + /// Resuelve ΔP(q) = dpTarget por Newton 1D con amortiguación. + /// + public double SolveFlowGivenDp(double dpTarget, Fluid fluid) + { + var q = Q; + var qMax = 1.0; // m³/s máximo razonable + + for (int i = 0; i < 50; i++) // más iteraciones + { + var f = Dp(q, fluid) - dpTarget; + var df = DdpDq(q, fluid); + + // Verificar división por cero + if (Math.Abs(df) < 1e-20) + df = 1e-10; + + var step = f / df; + + // Limitar el paso para estabilidad + var maxStep = Math.Min(0.1, Math.Abs(q) * 0.5 + 0.01); + if (Math.Abs(step) > maxStep) + step = Math.Sign(step) * maxStep; + + var qNew = q - step; + + // Limitar caudal dentro de rango razonable + qNew = Math.Max(-qMax, Math.Min(qMax, qNew)); + + // Amortiguación más conservadora + var relax = i < 10 ? 0.3 : 0.1; + q = (1 - relax) * q + relax * qNew; + + // Criterio de convergencia más estricto + if (Math.Abs(f) < 1e-1) // Pa + break; + + // Debug: evitar divergencia extrema + if (Math.Abs(q) > qMax) + q = Math.Sign(q) * qMax * 0.1; + } + + Q = q; + return q; + } + } +} diff --git a/HydraulicSimulator/Models/DischargeTank.cs b/HydraulicSimulator/Models/DischargeTank.cs new file mode 100644 index 0000000..2845eaa --- /dev/null +++ b/HydraulicSimulator/Models/DischargeTank.cs @@ -0,0 +1,79 @@ +using System; +using System.Linq; +using System.Collections.Generic; + +namespace HydraulicSimulator.Models +{ + /// + /// Tanque especial para descarga libre con cálculo de nivel dinámico + /// + public class DischargeTank : Element + { + public double Area { get; set; } // m² - área del tanque + public double CurrentVolume { get; set; } // m³ - volumen actual + public double MaxVolume { get; set; } // m³ - volumen máximo + public double MinVolume { get; set; } // m³ - volumen mínimo + public string TankId { get; set; } // ID del tanque + + public DischargeTank(string tankId, double area, double initialVolume = 0, + double maxVolume = double.MaxValue, double minVolume = 0) + { + TankId = tankId; + Area = area; + CurrentVolume = initialVolume; + MaxVolume = maxVolume; + MinVolume = minVolume; + } + + /// + /// Altura actual del líquido en el tanque + /// + public double CurrentHeight => CurrentVolume / Area; + + /// + /// Presión hidrostática en el fondo del tanque + /// + public double BottomPressure(Fluid fluid) => fluid.Rho * 9.81 * CurrentHeight; + + /// + /// Actualizar volumen basado en flujo neto + /// + public void UpdateVolume(double netFlowRate, double deltaTime) + { + double volumeChange = netFlowRate * deltaTime; + CurrentVolume += volumeChange; + + // Limitar entre mínimo y máximo + CurrentVolume = Math.Max(MinVolume, Math.Min(MaxVolume, CurrentVolume)); + } + + /// + /// Para tanque de descarga, la caída de presión es mínima + /// + public override double Dp(double q, Fluid fluid) + { + // Caída de presión mínima para tanque de descarga + return Math.Sign(q) * 100.0; // 100 Pa de pérdida nominal + } + + public override double DdpDq(double q, Fluid fluid) + { + return 1e-6; // Derivada muy pequeña para estabilidad + } + + /// + /// Verificar si el tanque está desbordando + /// + public bool IsOverflowing => CurrentVolume >= MaxVolume; + + /// + /// Verificar si el tanque está vacío + /// + public bool IsEmpty => CurrentVolume <= MinVolume; + + /// + /// Porcentaje de llenado (0.0 - 1.0) + /// + public double FillPercentage => (CurrentVolume - MinVolume) / (MaxVolume - MinVolume); + } +} diff --git a/HydraulicSimulator/Models/Element.cs b/HydraulicSimulator/Models/Element.cs new file mode 100644 index 0000000..dc261b0 --- /dev/null +++ b/HydraulicSimulator/Models/Element.cs @@ -0,0 +1,21 @@ +using System; + +namespace HydraulicSimulator.Models +{ + /// + /// Clase base para todos los elementos hidráulicos + /// + public abstract class Element + { + /// + /// Delta de presión (Pa) con signo positivo en sentido de q (opone el flujo). + /// Las bombas devuelven negativo (agregan presión). + /// + public abstract double Dp(double q, Fluid fluid); + + /// + /// Derivada d(ΔP)/dq, usada en Newton 1D de rama. + /// + public abstract double DdpDq(double q, Fluid fluid); + } +} diff --git a/HydraulicSimulator/Models/Fluid.cs b/HydraulicSimulator/Models/Fluid.cs new file mode 100644 index 0000000..c43f7fe --- /dev/null +++ b/HydraulicSimulator/Models/Fluid.cs @@ -0,0 +1,41 @@ +using System; + +namespace HydraulicSimulator.Models +{ + /// + /// Representa un fluido con sus propiedades físicas + /// + public class Fluid + { + public double Rho { get; set; } = 998.0; // kg/m³ - densidad + public double Mu { get; set; } = 1.0e-3; // Pa*s - viscosidad dinámica + public string Name { get; set; } = "Water"; // nombre del fluido + public double Temp { get; set; } = 20.0; // °C - temperatura + public double PVapor { get; set; } = 2337.0; // Pa - presión de vapor a 20°C + + public Fluid(double rho = 998.0, double mu = 1.0e-3, string name = "Water", + double temp = 20.0, double pVapor = 2337.0) + { + Rho = rho; + Mu = mu; + Name = name; + Temp = temp; + PVapor = pVapor; + } + + /// + /// Viscosidad cinemática (m²/s) + /// + public double Nu => Mu / Rho; + + /// + /// Verifica si hay cavitación + /// + public bool IsCavitating(double pAbs) => pAbs < PVapor; + + // Fluidos predefinidos + public static Fluid Water20C => new Fluid(998.0, 1.0e-3, "Water_20C", 20.0); + public static Fluid Water60C => new Fluid(983.0, 4.66e-4, "Water_60C", 60.0); + public static Fluid OilSAE30 => new Fluid(876.0, 0.1, "Oil_SAE30", 20.0); + } +} diff --git a/HydraulicSimulator/Models/HydraulicNetwork.cs b/HydraulicSimulator/Models/HydraulicNetwork.cs new file mode 100644 index 0000000..c82a8ce --- /dev/null +++ b/HydraulicSimulator/Models/HydraulicNetwork.cs @@ -0,0 +1,312 @@ +using System; +using System.Collections.Generic; +using System.Linq; + +namespace HydraulicSimulator.Models +{ + /// + /// Resultado de la simulación de la red hidráulica + /// + public class SolutionResult + { + public bool Converged { get; set; } + public Dictionary Flows { get; set; } = new Dictionary(); + public Dictionary Pressures { get; set; } = new Dictionary(); + public int Iterations { get; set; } + public double Residual { get; set; } + public string ErrorMessage { get; set; } = ""; + + public SolutionResult(bool converged = false) + { + Converged = converged; + } + } + + /// + /// Red hidráulica y solver + /// + public class HydraulicNetwork + { + public Fluid Fluid { get; set; } + public Dictionary Nodes { get; set; } = new Dictionary(); + public List Branches { get; set; } = new List(); + + public HydraulicNetwork(Fluid fluid = null) + { + Fluid = fluid ?? Fluid.Water20C; + } + + public void AddNode(string name, double? fixedP = null) + { + if (Nodes.ContainsKey(name)) + return; + + if (fixedP == null) + { + Nodes[name] = new Node(name, fixedP: false, p: 0.0); + } + else + { + Nodes[name] = new Node(name, fixedP: true, p: fixedP.Value); + } + } + + public void AddElement(Element element, string n1, string n2, string name = "") + { + if (!Nodes.ContainsKey(n1) || !Nodes.ContainsKey(n2)) + throw new ArgumentException("Nodos inexistentes"); + + var branchName = string.IsNullOrEmpty(name) ? $"{n1}->{n2}" : name; + Branches.Add(new Branch(n1, n2, new List { element }, branchName)); + } + + public void AddBranch(string n1, string n2, List elements, string name = "") + { + if (!Nodes.ContainsKey(n1) || !Nodes.ContainsKey(n2)) + throw new ArgumentException("Nodos inexistentes"); + + Branches.Add(new Branch(n1, n2, elements, name)); + } + + public SolutionResult Solve(int maxIterations = 100, double tolerance = 1e-3, + double relaxationFactor = 0.1, bool verbose = false) + { + try + { + var names = Nodes.Keys.ToList(); + var free = names.Where(n => !Nodes[n].FixedP).ToList(); + var fixedNodes = names.Where(n => Nodes[n].FixedP).ToList(); + var idxFree = free.Select((n, i) => new { n, i }).ToDictionary(x => x.n, x => x.i); + + // Inicial: presiones ~ promedio de fijos + double pRef = 0.0; + if (fixedNodes.Any()) + { + pRef = fixedNodes.Average(n => Nodes[n].P); + } + + foreach (var n in free) + { + Nodes[n].P = pRef; + } + + // Inicializar caudales a valores pequeños + foreach (var b in Branches) + { + b.Q = 0.001; // m³/s inicial pequeño + } + + double normR = 0.0; + int it = 0; + + // Iteración global sobre presiones nodales + for (it = 0; it < maxIterations; it++) + { + // 1) con presiones actuales, resolvés q de cada rama + foreach (var b in Branches) + { + var dpTarget = Nodes[b.N1].P - Nodes[b.N2].P; + b.SolveFlowGivenDp(dpTarget, Fluid); + } + + // 2) residuos de continuidad en nodos libres + var R = new double[free.Count]; + var G = new double[free.Count, free.Count]; // Jacobiano nodal + + foreach (var b in Branches) + { + var i = b.N1; + var j = b.N2; + // sensibilidad dq/d(dp) = 1 / d(ΔP)/dq + var ddpDq = b.DdpDq(b.Q, Fluid); + + // Evitar división por cero y valores demasiado grandes + ddpDq = Math.Max(1e-10, Math.Min(1e10, Math.Abs(ddpDq))); + var k = 1.0 / ddpDq; + + // aporte a residuos (signo: positiva saliendo de nodo n1) + if (idxFree.ContainsKey(i)) + { + var idx_i = idxFree[i]; + R[idx_i] += b.Q; + G[idx_i, idx_i] += k; + if (idxFree.ContainsKey(j)) + { + var idx_j = idxFree[j]; + G[idx_i, idx_j] -= k; + } + } + if (idxFree.ContainsKey(j)) + { + var idx_j = idxFree[j]; + R[idx_j] -= b.Q; + G[idx_j, idx_j] += k; + if (idxFree.ContainsKey(i)) + { + var idx_i = idxFree[i]; + G[idx_j, idx_i] -= k; + } + } + } + + normR = R.Length > 0 ? R.Max(Math.Abs) : 0.0; + if (verbose) + Console.WriteLine($"it {it}: |R|_inf={normR:E3}"); + + if (normR < tolerance) + break; + + if (free.Count > 0) + { + // 3) resolver actualización de presiones + var dpUpdate = SolveLinearSystem(G, R, free.Count); + + // Limitar la actualización de presión + var maxDp = 50000.0; // 50 kPa máximo por iteración + for (int k = 0; k < dpUpdate.Length; k++) + { + dpUpdate[k] = Math.Max(-maxDp, Math.Min(maxDp, dpUpdate[k])); + } + + // amortiguación más conservadora + var relax = it < 10 ? relaxationFactor : relaxationFactor * 0.5; + for (int k = 0; k < dpUpdate.Length; k++) + { + dpUpdate[k] *= relax; + } + + foreach (var kvp in idxFree) + { + Nodes[kvp.Key].P += dpUpdate[kvp.Value]; + } + } + } + + // Preparar resultado + var result = new SolutionResult(normR < tolerance) + { + Iterations = it, + Residual = normR + }; + + // Llenar flows y pressures + foreach (var b in Branches) + { + result.Flows[b.Name] = b.Q; + } + + foreach (var kvp in Nodes) + { + result.Pressures[kvp.Key] = kvp.Value.P; + } + + return result; + } + catch (Exception ex) + { + return new SolutionResult(false) + { + ErrorMessage = ex.Message, + Iterations = 0, + Residual = double.MaxValue + }; + } + } + + private double[] SolveLinearSystem(double[,] G, double[] R, int n) + { + // Implementación simple de eliminación gaussiana con pivoteo parcial + var A = new double[n, n + 1]; + + // Copiar G y -R al sistema aumentado + for (int i = 0; i < n; i++) + { + for (int j = 0; j < n; j++) + { + A[i, j] = G[i, j]; + } + A[i, n] = -R[i]; + } + + // Regularización + for (int i = 0; i < n; i++) + { + A[i, i] += 1e-6; + } + + // Eliminación gaussiana + for (int i = 0; i < n; i++) + { + // Encontrar el pivote + int maxRow = i; + for (int k = i + 1; k < n; k++) + { + if (Math.Abs(A[k, i]) > Math.Abs(A[maxRow, i])) + maxRow = k; + } + + // Intercambiar filas + for (int k = i; k <= n; k++) + { + var temp = A[maxRow, k]; + A[maxRow, k] = A[i, k]; + A[i, k] = temp; + } + + // Hacer todos los elementos debajo del pivote igual a 0 + for (int k = i + 1; k < n; k++) + { + if (Math.Abs(A[i, i]) < 1e-12) + A[i, i] = 1e-12; + + var c = A[k, i] / A[i, i]; + for (int j = i; j <= n; j++) + { + if (i == j) + A[k, j] = 0; + else + A[k, j] -= c * A[i, j]; + } + } + } + + // Sustitución hacia atrás + var solution = new double[n]; + for (int i = n - 1; i >= 0; i--) + { + solution[i] = A[i, n]; + for (int j = i + 1; j < n; j++) + { + solution[i] -= A[i, j] * solution[j]; + } + + if (Math.Abs(A[i, i]) < 1e-12) + A[i, i] = 1e-12; + + solution[i] = solution[i] / A[i, i]; + } + + return solution; + } + + /// + /// Genera un reporte de la solución + /// + public void Report() + { + Console.WriteLine("== Nodos (Pa) =="); + foreach (var kvp in Nodes) + { + var node = kvp.Value; + var kind = node.FixedP ? "FIX" : "FREE"; + Console.WriteLine($"{node.Name,10}: {node.P,12:F1} [{kind}]"); + } + + Console.WriteLine("\n== Ramas (q m³/s) =="); + foreach (var b in Branches) + { + Console.WriteLine($"{b.Name,15}: {b.Q,10:E6}"); + } + } + } +} diff --git a/HydraulicSimulator/Models/MinorLoss.cs b/HydraulicSimulator/Models/MinorLoss.cs new file mode 100644 index 0000000..dae484a --- /dev/null +++ b/HydraulicSimulator/Models/MinorLoss.cs @@ -0,0 +1,35 @@ +using System; + +namespace HydraulicSimulator.Models +{ + /// + /// Pérdida menor con coeficiente K + /// + public class MinorLoss : Element + { + public double K { get; set; } // adimensional + public double DRef { get; set; } // m (para el área de referencia) + + public MinorLoss(double k, double dRef) + { + K = k; + DRef = dRef; + } + + private double Area => Math.PI * (DRef * DRef) / 4.0; + + public override double Dp(double q, Fluid fluid) + { + var area = Area; + var c = K * 0.5 * fluid.Rho / (area * area); + return c * q * Math.Abs(q); + } + + public override double DdpDq(double q, Fluid fluid) + { + var area = Area; + var c = K * 0.5 * fluid.Rho / (area * area); + return 2.0 * c * Math.Abs(q) + 1e-12; + } + } +} diff --git a/HydraulicSimulator/Models/Node.cs b/HydraulicSimulator/Models/Node.cs new file mode 100644 index 0000000..f1989ca --- /dev/null +++ b/HydraulicSimulator/Models/Node.cs @@ -0,0 +1,19 @@ +namespace HydraulicSimulator.Models +{ + /// + /// Representa un nodo en la red hidráulica + /// + public class Node + { + public string Name { get; set; } + public bool FixedP { get; set; } = false; + public double P { get; set; } = 0.0; // Pa + + public Node(string name, bool fixedP = false, double p = 0.0) + { + Name = name; + FixedP = fixedP; + P = p; + } + } +} diff --git a/HydraulicSimulator/Models/Pipe.cs b/HydraulicSimulator/Models/Pipe.cs new file mode 100644 index 0000000..0cda07f --- /dev/null +++ b/HydraulicSimulator/Models/Pipe.cs @@ -0,0 +1,53 @@ +using System; + +namespace HydraulicSimulator.Models +{ + /// + /// Tubería con ecuación de Darcy-Weisbach y factor de fricción por Swamee-Jain + /// + public class Pipe : Element + { + public double L { get; set; } // m - longitud + public double D { get; set; } // m - diámetro + public double Rough { get; set; } = 4.5e-5; // m - rugosidad (acero comercial ~45 micrones) + + public Pipe(double length, double diameter, double roughness = 4.5e-5) + { + L = length; + D = diameter; + Rough = roughness; + } + + private double Area => Math.PI * (D * D) / 4.0; + + private double FrictionFactor(double q, Fluid fluid) + { + var area = Area; + var v = Math.Abs(q) / area; + var re = Math.Max(4000.0, fluid.Rho * v * D / fluid.Mu); // forzamos turbulento + var epsRel = Rough / D; + // Swamee-Jain + var f = 0.25 / Math.Pow(Math.Log10(epsRel / 3.7 + 5.74 / Math.Pow(re, 0.9)), 2); + return f; + } + + public override double Dp(double q, Fluid fluid) + { + var area = Area; + var f = FrictionFactor(q, fluid); + var k = f * (L / D) * 0.5 * fluid.Rho / (area * area); + return k * q * Math.Abs(q); // signo + } + + public override double DdpDq(double q, Fluid fluid) + { + // Ignoramos df/dq para estabilidad/simplicidad (funciona muy bien). + var area = Area; + var qAbs = Math.Max(Math.Abs(q), 1e-9); + var qSign = q >= 0 ? 1 : -1; + var f = FrictionFactor(qAbs * qSign, fluid); + var k = f * (L / D) * 0.5 * fluid.Rho / (area * area); + return 2.0 * k * Math.Abs(q) + 1e-12; // evita 0 + } + } +} diff --git a/HydraulicSimulator/Models/PumpHQ.cs b/HydraulicSimulator/Models/PumpHQ.cs new file mode 100644 index 0000000..cb782d5 --- /dev/null +++ b/HydraulicSimulator/Models/PumpHQ.cs @@ -0,0 +1,50 @@ +using System; + +namespace HydraulicSimulator.Models +{ + /// + /// Bomba con curva H(Q)=H0*(1-(Q/Q0)²) y ley de afinidad con velocidad relativa + /// + public class PumpHQ : Element + { + public double H0 { get; set; } // m, a velocidad nominal (shutoff head) + public double Q0 { get; set; } // m³/s, caudal a cabeza cero, vel nominal + public double SpeedRel { get; set; } = 1.0; // n / n_nominal + public int Direction { get; set; } = 1; // +1 si impulsa de i->j, -1 si al revés + + public PumpHQ(double h0, double q0, double speedRel = 1.0, int direction = 1) + { + H0 = h0; + Q0 = q0; + SpeedRel = speedRel; + Direction = direction; + } + + private (double H0s, double Q0s) Scaled + { + get + { + var s = Math.Max(1e-3, SpeedRel); + return (H0 * (s * s), Q0 * s); + } + } + + public override double Dp(double q, Fluid fluid) + { + var (h0s, q0s) = Scaled; + // Limitamos fuera de rango para estabilidad + var qq = Math.Max(-q0s * 0.999, Math.Min(q0s * 0.999, q)); + var h = h0s * (1.0 - Math.Pow(qq / q0s, 2)); + var dp = -Direction * fluid.Rho * 9.80665 * h; + // dp es negativo si la bomba agrega presión en el sentido de la rama + return dp; + } + + public override double DdpDq(double q, Fluid fluid) + { + var (h0s, q0s) = Scaled; + var dhDq = -2.0 * h0s * q / (q0s * q0s); + return -Direction * fluid.Rho * 9.80665 * dhDq + 1e-12; + } + } +} diff --git a/HydraulicSimulator/Models/ValveKv.cs b/HydraulicSimulator/Models/ValveKv.cs new file mode 100644 index 0000000..51356af --- /dev/null +++ b/HydraulicSimulator/Models/ValveKv.cs @@ -0,0 +1,44 @@ +using System; + +namespace HydraulicSimulator.Models +{ + /// + /// Válvula por Kv con apertura 0..1 + /// + public class ValveKv : Element + { + public double KvFull { get; set; } // m³/h / sqrt(bar) a 100% + public double Opening { get; set; } // 0..1 + + public ValveKv(double kvFull, double opening) + { + KvFull = kvFull; + Opening = opening; + } + + private double KvEff + { + get + { + var x = Math.Min(1.0, Math.Max(0.0, Opening)); + return Math.Max(1e-6, KvFull * x); // lineal simple y evita 0 + } + } + + public override double Dp(double q, Fluid fluid) + { + var kv = KvEff; + var qh = q * 3600.0; // m³/h + var dpBar = Math.Pow(qh / kv, 2); // bar (SG≈1) + var dpPa = dpBar * 1e5; + return Math.Sign(q) * Math.Abs(dpPa); + } + + public override double DdpDq(double q, Fluid fluid) + { + var kv = KvEff; + var coeff = 2.0 * 1e5 * Math.Pow(3600.0 / kv, 2); + return coeff * Math.Abs(q) + 1e-12; + } + } +}