Base antes de agregar logica de HydraulicSimulator

This commit is contained in:
Miguel 2025-09-04 11:41:26 +02:00
parent 4f2a109332
commit 3fe5b5497f
12 changed files with 1134 additions and 0 deletions

View File

@ -188,6 +188,7 @@
</ItemGroup>
<ItemGroup>
<Folder Include="ObjetosSim\HydraulicSimulator\" />
<Folder Include="paddleocr\cls\inference\" />
<Folder Include="paddleocr\det\inference\" />
<Folder Include="paddleocr\keys\" />

View File

@ -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<Element> { pump, mainPipe });
network.AddBranch("UNION_1", "PROCESO_1",
new List<Element> { valve });
network.AddBranch("PROCESO_1", "DESCARGA",
new List<Element> { 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<GraphicElement> 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<double>(graphic, "Head"),
q0: GetProperty<double>(graphic, "MaxFlow") / 3600.0 // Convertir m³/h a m³/s
),
"Pipe" => new Pipe(
length: GetProperty<double>(graphic, "Length"),
diameter: GetProperty<double>(graphic, "Diameter"),
roughness: GetProperty<double>(graphic, "Roughness", 0.0015)
),
"Valve" => new ValveKv(
kvFull: GetProperty<double>(graphic, "Kv"),
opening: GetProperty<double>(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> { element },
graphic.Id
);
}
return network;
}
private T GetProperty<T>(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<string, double> 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<Element> { pump1, supply1, valve1 });
network.AddBranch("TANQUE_PRINCIPAL", "MIXER_2",
new List<Element> { pump2, supply2, valve2 });
network.AddBranch("MIXER_1", "DESCARGA",
new List<Element> { discharge1 });
network.AddBranch("MIXER_2", "DESCARGA",
new List<Element> { 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?

View File

@ -0,0 +1,82 @@
using System;
using System.Collections.Generic;
using System.Linq;
namespace HydraulicSimulator.Models
{
/// <summary>
/// Rama (elementos en serie)
/// </summary>
public class Branch
{
public string N1 { get; set; }
public string N2 { get; set; }
public List<Element> 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<Element> elements, string name = "")
{
N1 = n1;
N2 = n2;
Elements = new List<Element>(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));
}
/// <summary>
/// Resuelve ΔP(q) = dpTarget por Newton 1D con amortiguación.
/// </summary>
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;
}
}
}

View File

@ -0,0 +1,79 @@
using System;
using System.Linq;
using System.Collections.Generic;
namespace HydraulicSimulator.Models
{
/// <summary>
/// Tanque especial para descarga libre con cálculo de nivel dinámico
/// </summary>
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;
}
/// <summary>
/// Altura actual del líquido en el tanque
/// </summary>
public double CurrentHeight => CurrentVolume / Area;
/// <summary>
/// Presión hidrostática en el fondo del tanque
/// </summary>
public double BottomPressure(Fluid fluid) => fluid.Rho * 9.81 * CurrentHeight;
/// <summary>
/// Actualizar volumen basado en flujo neto
/// </summary>
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));
}
/// <summary>
/// Para tanque de descarga, la caída de presión es mínima
/// </summary>
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
}
/// <summary>
/// Verificar si el tanque está desbordando
/// </summary>
public bool IsOverflowing => CurrentVolume >= MaxVolume;
/// <summary>
/// Verificar si el tanque está vacío
/// </summary>
public bool IsEmpty => CurrentVolume <= MinVolume;
/// <summary>
/// Porcentaje de llenado (0.0 - 1.0)
/// </summary>
public double FillPercentage => (CurrentVolume - MinVolume) / (MaxVolume - MinVolume);
}
}

View File

@ -0,0 +1,21 @@
using System;
namespace HydraulicSimulator.Models
{
/// <summary>
/// Clase base para todos los elementos hidráulicos
/// </summary>
public abstract class Element
{
/// <summary>
/// Delta de presión (Pa) con signo positivo en sentido de q (opone el flujo).
/// Las bombas devuelven negativo (agregan presión).
/// </summary>
public abstract double Dp(double q, Fluid fluid);
/// <summary>
/// Derivada d(ΔP)/dq, usada en Newton 1D de rama.
/// </summary>
public abstract double DdpDq(double q, Fluid fluid);
}
}

View File

@ -0,0 +1,41 @@
using System;
namespace HydraulicSimulator.Models
{
/// <summary>
/// Representa un fluido con sus propiedades físicas
/// </summary>
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;
}
/// <summary>
/// Viscosidad cinemática (m²/s)
/// </summary>
public double Nu => Mu / Rho;
/// <summary>
/// Verifica si hay cavitación
/// </summary>
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);
}
}

View File

@ -0,0 +1,312 @@
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));
}
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;
}
/// <summary>
/// Genera un reporte de la solución
/// </summary>
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}");
}
}
}
}

View File

@ -0,0 +1,35 @@
using System;
namespace HydraulicSimulator.Models
{
/// <summary>
/// Pérdida menor con coeficiente K
/// </summary>
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;
}
}
}

View File

@ -0,0 +1,19 @@
namespace HydraulicSimulator.Models
{
/// <summary>
/// Representa un nodo en la red hidráulica
/// </summary>
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;
}
}
}

View File

@ -0,0 +1,53 @@
using System;
namespace HydraulicSimulator.Models
{
/// <summary>
/// Tubería con ecuación de Darcy-Weisbach y factor de fricción por Swamee-Jain
/// </summary>
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
}
}
}

View File

@ -0,0 +1,50 @@
using System;
namespace HydraulicSimulator.Models
{
/// <summary>
/// Bomba con curva H(Q)=H0*(1-(Q/Q0)²) y ley de afinidad con velocidad relativa
/// </summary>
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;
}
}
}

View File

@ -0,0 +1,44 @@
using System;
namespace HydraulicSimulator.Models
{
/// <summary>
/// Válvula por Kv con apertura 0..1
/// </summary>
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;
}
}
}