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)); } /// /// Actualiza las presiones en bombas con verificación de NPSH /// private void UpdatePumpPressures() { var currentPressures = new Dictionary(); foreach (var kvp in Nodes) { currentPressures[kvp.Key] = kvp.Value.P; } foreach (var branch in Branches) { foreach (var element in branch.Elements) { if (element is PumpHQWithSuctionCheck pumpWithCheck) { pumpWithCheck.UpdatePressures(currentPressures); } } } } 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++) { // Actualizar presiones en bombas con verificación de NPSH UpdatePumpPressures(); // 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; // Console output deshabilitado para mejorar rendimiento // 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() { // Reporte deshabilitado para mejorar rendimiento /* 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}"); } */ } } }