343 lines
12 KiB
C#
343 lines
12 KiB
C#
using System;
|
|
using System.Collections.Generic;
|
|
using System.Linq;
|
|
|
|
namespace HydraulicSimulator.Models
|
|
{
|
|
/// <summary>
|
|
/// Resultado de la simulación de la red hidráulica
|
|
/// </summary>
|
|
public class SolutionResult
|
|
{
|
|
public bool Converged { get; set; }
|
|
public Dictionary<string, double> Flows { get; set; } = new Dictionary<string, double>();
|
|
public Dictionary<string, double> Pressures { get; set; } = new Dictionary<string, double>();
|
|
public int Iterations { get; set; }
|
|
public double Residual { get; set; }
|
|
public string ErrorMessage { get; set; } = "";
|
|
|
|
public SolutionResult(bool converged = false)
|
|
{
|
|
Converged = converged;
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Red hidráulica y solver
|
|
/// </summary>
|
|
public class HydraulicNetwork
|
|
{
|
|
public Fluid Fluid { get; set; }
|
|
public Dictionary<string, Node> Nodes { get; set; } = new Dictionary<string, Node>();
|
|
public List<Branch> Branches { get; set; } = new List<Branch>();
|
|
|
|
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> { element }, branchName));
|
|
}
|
|
|
|
public void AddBranch(string n1, string n2, List<Element> elements, string name = "")
|
|
{
|
|
if (!Nodes.ContainsKey(n1) || !Nodes.ContainsKey(n2))
|
|
throw new ArgumentException("Nodos inexistentes");
|
|
|
|
Branches.Add(new Branch(n1, n2, elements, name));
|
|
}
|
|
|
|
/// <summary>
|
|
/// Actualiza las presiones en bombas con verificación de NPSH
|
|
/// </summary>
|
|
private void UpdatePumpPressures()
|
|
{
|
|
var currentPressures = new Dictionary<string, double>();
|
|
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;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Genera un reporte de la solución
|
|
/// </summary>
|
|
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}");
|
|
}
|
|
*/
|
|
}
|
|
}
|
|
}
|