tmabit
/
BUDEES
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Initial commit

This commit is contained in:
Thomas Mabit 2024-11-19 12:50:18 +01:00
commit bc81ddd74e
21 changed files with 3127 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
src/Main.cpp Normal file
View File

@ -0,0 +1,56 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
struct JNICaller {};
#include "Simulateur/Simulateur.h"
int main() {
// -------- Parameters -------- //
Configuration configuration;
// System dimensions [m]
configuration.L[0] = 2.64e-3;
configuration.L[1] = 4.2e-3;
// Current density [A/m2]
configuration.J = 500;
// Simulation duration [s]
configuration.tf = 4.0;//24.25;
// --------- Topology --------- //
Topology topology;
// Polygon count
topology.polygonCount[0] = 1;
topology.polygonCount[1] = 1;
// Polygon indexs
topology.polygonIndex[0] = new int[] {0, 6};
topology.polygonIndex[1] = new int[] {6, 12};
// Polygon defs
Point points[] = {
{0.0, 0.0}, {0.27, 0.0}, {0.44, 0.63}, {0.43, 0.8}, {0.27, 1.0}, {0.0, 1.0},
{1.0, 0.0}, {1.0, 1.0}, {0.60, 1.0}, {0.72, 0.8}, {0.71, 0.46}, {0.60, 0.0},
}; topology.polygonDefinitions = points;
// Simulation
if (!fs::exists(outputFolder)) {fs::create_directory(outputFolder);} outputDelete(outputFolder.c_str());
Performances performances = simulate(configuration, topology, nullptr);
cout << endl << "Simulation complete (" << (performances.collected[H2]+performances.lost[H2]+performances.remaining[H2]) << "mol of H2 produced) (n=" << (performances.collected[H2] / (performances.collected[H2]+performances.lost[H2])) << "%, R=" << performances.resistor << "ohm.m2, p=" << (performances.collected[H2] / (performances.collected[H2]+performances.contamination[H2])) << ") !" << endl;
// Cleanup memory
delete[] topology.polygonIndex[0];
delete[] topology.polygonIndex[1];
// Exit
return 0;
}

157
src/Simulateur/Simulateur.h Normal file
View File

@ -0,0 +1,157 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
#include <iostream>
#include <stdlib.h>
#include <vector>
#include <list>
#include <chrono>
#include <fstream>
#include <string>
#include <filesystem>
#include <stdexcept>
using namespace std;
#include <loadbmp/loadbmp.h>
#include <Eigen/Dense>
#include <Eigen/SparseLU>
#include <Eigen/SparseCholesky>
typedef Eigen::MatrixXd Mat;
typedef Eigen::SparseMatrix<double> SpMat;
typedef Eigen::Triplet<double> Trip;
using Eigen::VectorXd;
using Eigen::MatrixXd;
using Eigen::MatrixXi;
#include "params/params.h"
#include "utils/typesdef.h"
#include "utils/math2.h"
#include "params/consts.h"
#include "params/loader.h"
#include "utils/polygon.h"
#include "utils/progress.h"
#include "solver/Tension.h"
#include "solver/Current.h"
#include "solver/Topology.h"
#include "solver/Nucleation.h"
#include "solver/Timestep.h"
#include "solver/Dynamics.h"
#include "solver/Counter.h"
#include "utils/output.h"
#include "utils/loader.h"
void capture(Configuration& configuration, Topology& topology, const char *path) {
// Load parameters and output mask
Parameters parameters = loadParameters(configuration);
MatrixXi MASK = generateMask(parameters, topology);
outputMatrixI(path, MASK, 0, 2);
}
Performances simulate(Configuration& configuration, Topology& topology, JNICaller* jniCaller) {
// Load parameters and compute pre-valus
Parameters parameters = loadParameters(configuration);
// Prepare progressbar
progressInit(jniCaller);
// Electrical simulation
MatrixXi MASK = generateMask(parameters, topology);
MatrixTension TENSION = solveTension(parameters, MASK, false);
MatrixXd CURRENT = computeCurrent(parameters, MASK, TENSION);
MatrixXd SOURCES = computeSourceNormalized(parameters, MASK, TENSION);
vector<Site> Sites = buildNucleationSite(parameters, topology, SOURCES);
// Bubble dynamics initialization
TimeSteps timesteps = {0, 0, 0, 0};
ListBubbles Bubbles; Bubbles.flagAttached = false; Bubbles.flagDetached = false;
Bubbles.countAttached = 0; Bubbles.countDetached = 0;
Bubbles.tNextNucleation = 0;
vector<Wall> Walls = generateWalls(parameters, topology);
ElectrolysisStatistics Statistics = {{0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, 0.0, 0.0, {0.0, 0.0}, {0.0, 0.0}};
updateStatistics(Statistics, TENSION, CURRENT);
// Output electrical (only if not called from jni)
if (outputElectrical && (jniCaller == nullptr)) {
outputMatrixI(outputFolder, "/mask.bmp", MASK, -1, 1);
outputMatrix(outputFolder, "/tension.bmp", TENSION.VALUES, TENSION.min, TENSION.max, -1.0);
outputData(outputFolder, "/tension.csv", parameters, TENSION.VALUES, TENSION.min, TENSION.max);
outputMatrix(outputFolder, "/current.bmp", CURRENT, -CURRENT.maxCoeff(), CURRENT.maxCoeff(), -1.0);
outputPoints(outputFolder, "/source.bmp", SOURCES, MAX(-SOURCES.minCoeff(), SOURCES.maxCoeff()), 7);
outputSites(outputFolder, "/sites.bmp", parameters, Sites);
}
// Bubble dynamics simulation
double t = 0.0; bool stopFlag = false;
while ((t < parameters.tf) && (!stopFlag)) {
// 1| process nucleation
nucleateBubble(Bubbles, Sites, t);
// 2| process bubble evolution
computeBubble(parameters, Bubbles, Walls, Statistics, t, timesteps.dt, &stopFlag);
// 3| update tension
if (electricalCoupling) if (t > TENSION.tNext) {
MatrixXi MASK_new = updateMask(parameters, MASK, Bubbles);
if (MASK != MASK_new) {
MASK = MASK_new;
MatrixTension TENSION_new = solveTension(parameters, MASK, true);
if (TENSION_new.sucess) {
TENSION = TENSION_new;
CURRENT = computeCurrent(parameters, MASK, TENSION);
SOURCES = computeSourceNormalized(parameters, MASK, TENSION);
updateNucleationSite(parameters, SOURCES, Sites);
updateStatistics(Statistics, TENSION, CURRENT);
}
TENSION.tNext += timesteps.dtElec;
}
}
// 4| update timestep
t += timesteps.dt; timesteps.i++;
timesteps = updateTimestep(timesteps, Bubbles.maxRadius, Bubbles.maxVelocity, Sites, &stopFlag);
// 5| output
progressRefresh1((t/parameters.tf), jniCaller);
if (jniCaller == nullptr) {
// Frame output
outputFrame(outputFolder, "/frame.blb", t, parameters, timesteps, topology, Bubbles, Statistics);
// Electrical outputs
if (outputElectrical && electricalCoupling) {
timesteps.iElec++;
outputMatrixI(outputFolder, "/mask/", timesteps.iElec, MASK, -1, 1);
outputMatrix(outputFolder, "/tension/", timesteps.iElec, TENSION.VALUES, TENSION.min, TENSION.max, -1.0);
outputMatrix(outputFolder, "/current/", timesteps.iElec, CURRENT, -CURRENT.maxCoeff(), CURRENT.maxCoeff(), -1.0);
outputPoints(outputFolder, "/source/", timesteps.iElec, SOURCES, MAX(-SOURCES.minCoeff(), SOURCES.maxCoeff()), 7);
outputSites(outputFolder, "/sites/", timesteps.iElec, parameters, Sites);
}
}
}
// Performances
Performances performances;
if (!stopFlag) {
// Compute performances
performances = computePerformances(parameters, Statistics, Bubbles);
} else {
// Unable to compute performances
performances = zeroPerformances();
}
return performances;
}

92
src/Simulateur/params/consts.h Normal file
View File

@ -0,0 +1,92 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
/* Global const */
double Fcst = 96485.332; // Faraday cst [A.s/mol]
double Rcst = 8.314463; // Constante gaz parfait [J/mol.K]
/// Static pre-computation
double eThick, RSite[2], Fcapilaire[2], FcapilaireSq[2], RFritzExtra[2], Rmini[2],
nRPrefactor, nFritz[2], RFritzMax, RFritzMin, buoyancyPrefactor, QnPrefactor[2], Tprefactor[2], growingKineticFractionWait, uAscentMin,
tafelJ0A0[2], bPrefactor,
interactionLongMaximumRangeSq, interactionLubMaximumRangeSq, interactionRepulsiveRangeSq, interactionCollisionRangeSq,
interactionSafetyHalfRange, interactionSafetyHalfRangeInv, interactionSafetyHalfRangeSq, interactionSafetyRange, interactionSafetyRangeSq;
int bubbleLifeCount;
void staticPreComputation() {
// Diametre moyen
eThick = (RFritz[0] + RFritz[1]);
for (int i = 0; i < 2; i++) {
// Taille du défaut 'r' associé au site de nucléation pour le rayon de Fritz considérée (tel que Fcap = Fbuoyancy)
RSite[i] = 2*CUBE(RFritz[i])*rho*gravity/(3*gammaLiquidGas[i]);
// Force capilaire sur un site (Fcap = 2 PI r gamma)
Fcapilaire[i] = 2.0*M_PI*RSite[i]*gammaLiquidGas[i];
FcapilaireSq[i] = SQ(Fcapilaire[i]);
//...
RFritzExtra[i] = (extraRFritz+1)*RFritz[i];
Rmini[i] = detachementMinimalRadiusNorm*RFritz[i];
}
// Préfacteur pour le calcul de la quantité de matière dans une sphère (n = [4/3]PI * P0/RT * r^3)
nRPrefactor = 4.0*M_PI*P0/(3.0*Rcst*T0);
// Quantité de matière au détachement de Fritz
for (int i = 0; i < 2; i++) {
nFritz[i] = nRPrefactor*CUBE(RFritz[i]);
}
// ...
RFritzMax = mathMax(RFritz[0], RFritz[1]);
RFritzMin = mathMin(RFritz[0], RFritz[1]);
// Préfacteur pour la force de flottaison Fbuoyancy = [4/3]PI*rho*g * r^3
buoyancyPrefactor = (4.0/3.0)*M_PI*rho*gravity;
for (int i = 0; i < 2; i++) {
// Préfacteur pour le calcul de la charge (Jelec = [Nelec*F] * Jmolaire)
QnPrefactor[i] = Nelectronic[i]*Fcst;
// Préfacteur pour le calcul de la periode de nucleation
Tprefactor[i] = Nelectronic[i]*Fcst*nFritz[i];
}
// ...
growingKineticFractionWait = 1.0-growingKineticFraction;
// Plus petite vitesse d'ascension naturelle (pour les plus petites bulles c-a-d les bulles à RFritz)
uAscentMin = 2*rho*gravity*mathSQ(RFritzMin)/(9*mu);
// Préfacteur de la dérivée de Tafel
for (int i = 0; i < 2; i++) {
tafelJ0A0[i] = tafelJ0[i]/tafelA0[i];
}
// Préfacteur du terme de flotaison
bPrefactor = (2.0*rho*gravity)/(9.0*mu); // = (4/3)*PI*rho*g*R3 / 6*pi*mu*R
// Répétition des bulles -> 2 vies
bubbleLifeCount = (repeatYBubble ? 2 : 1);
// Version caré et inverse des parametres numériques
interactionLongMaximumRangeSq = mathSQ(interactionLongMaximumRange);
interactionLubMaximumRangeSq = mathSQ(interactionLubMaximumRange);
interactionRepulsiveRangeSq = mathSQ(interactionRepulsiveRange);
interactionCollisionRangeSq = mathSQ(interactionCollisionRange);
interactionSafetyHalfRange = interactionSafetyD+1.0;
interactionSafetyHalfRangeInv = 1.0/interactionSafetyHalfRange;
interactionSafetyHalfRangeSq = mathSQ(interactionSafetyHalfRange);
interactionSafetyRange = 2.0+interactionSafetyD;
interactionSafetyRangeSq = mathSQ(interactionSafetyRange);
}
#define TafelAnodeJ(deltaphi) (tafelJ0[ANODE]*exp((deltaphi)/tafelA0[ANODE]))
#define TafelCathodeJ(deltaphi) (-tafelJ0[CATHODE]*exp((deltaphi)/tafelA0[CATHODE]))
#define dTafelAnodeJ(deltaphi) (-tafelJ0A0[ANODE]*exp((deltaphi)/tafelA0[ANODE]))
#define dTafelCathodeJ(deltaphi) (-tafelJ0A0[CATHODE]*exp((deltaphi)/tafelA0[CATHODE]))

127
src/Simulateur/params/loader.h Normal file
View File

@ -0,0 +1,127 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
double str2bool(string value) {
return value == "true";
}
int str2int(string value) {
size_t pos;
int intValue = stoi(value, &pos);
if (pos != value.length()) throw invalid_argument("null");
return intValue;
}
double str2double(string value) {
size_t pos;
double doubleValue = stod(value, &pos);
if (pos != value.length()) throw invalid_argument("null");
return doubleValue;
}
void loadParamsFromFile(Configuration& configuration) {
ifstream file("params.txt");
if (file.is_open()) {
string line;
while (getline(file, line)) {
if (line.empty() || line[0] == '#') continue;
// Find '='
size_t pos = line.find("=");
if (pos == string::npos) {
cerr << "Error while parsing params.txt (" << line << ")" << endl << flush;
exit(1);
}
// Get variable and value
string variable = line.substr(0, pos), value = line.substr(pos + 1);
// Update
try {
if (variable == "dSite") dSite = str2double(value);
else if (variable == "RFritz") {RFritz[0] = str2double(value); RFritz[1] = str2double(value);}
else if (variable == "RFritz[0]") RFritz[0] = str2double(value);
else if (variable == "RFritz[1]") RFritz[1] = str2double(value);
else if (variable == "theta") {theta[0] = str2double(value); theta[1] = str2double(value);}
else if (variable == "theta[0]") theta[0] = str2double(value);
else if (variable == "theta[1]") theta[1] = str2double(value);
else if (variable == "sigma") sigma = str2double(value);
else if (variable == "electricalCoupling") electricalCoupling = str2bool(value);
else if (variable == "meshN") meshN = str2int(value);
else if (variable == "outputElectrical") outputElectrical = str2bool(value);
else if (variable == "outputFolder") outputFolder = value;
else if (variable == "useSubProgressBar") useSubProgressBar = str2bool(value);
else if (variable == "J") configuration.J = str2double(value);
else if (variable == "tf") configuration.tf = str2double(value);
else throw invalid_argument("null");
} catch (const invalid_argument& e) {
cerr << "Error while parsing params.txt (" << line << ")" << endl << flush;
exit(1);
}
}
}
}
// Dynamic pre-computation
struct Parameters {
double L[2]; /* = dynamic precomputation */ // [m]
int N[2]; /* = dynamic precomputation */
double dL[2]; /* = dynamic precomputation */ // [m]
double tf, tClogg; /* = dynamic precomputation */ // [s]
double J; /* = dynamic precomputation */ // [A/m2]
double I; /* = dynamic precomputation */ // [A]
};
Parameters loadParameters(Configuration& configuration) {
// Load params file and static precomputation
loadParamsFromFile(configuration);
staticPreComputation();
Parameters parameters;
// Dimensions
parameters.L[0] = configuration.L[0];
parameters.L[1] = configuration.L[1];
// Mesh
if (parameters.L[0] > parameters.L[1]) {
parameters.N[0] = meshN;
parameters.N[1] = round(((double) meshN)*(parameters.L[1]/parameters.L[0]));
} else {
parameters.N[1] = meshN;
parameters.N[0] = round(((double) meshN)*(parameters.L[0]/parameters.L[1]));
}
// Refine mesh
if (electricalCoupling) {
int refineFactor = 0;
for (int c = 0; c < 2; c++) {
int factor = (int) ceil(BUBBLE_MESH_COEFF * (parameters.L[c] / (2*Rmini[c])) / ((double) parameters.N[c]));
if (factor > 1) {
parameters.N[c] *= factor;
refineFactor = mathMax(refineFactor, factor);
}
}
}
// Compute dL/dR
parameters.dL[0] = configuration.L[0]/(parameters.N[0]-1);
parameters.dL[1] = configuration.L[1]/(parameters.N[1]-1);
// Current which correspond to the effective electrode surface (S = e x Ly)
parameters.J = configuration.J;
parameters.I = configuration.J*eThick*parameters.L[1];
// Others
parameters.tf = configuration.tf;
parameters.tClogg = cloggedTimeFactor * (parameters.L[1] / uAscentMin);
return parameters;
}

107
src/Simulateur/params/params.h Normal file
View File

@ -0,0 +1,107 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
#define ANODE 0
#define O2 0
#define CATHODE 1
#define H2 1
/* --------- Material --------- */
double dSite = 1.5/(1e-3); // Hole density [1/m]
double RFritz[2] = {60e-6, 60e-6}; // Rayon de Fritz [m]
double theta[2] = {M_PI, M_PI}; // Angle de mouillage des bulles (de l'ordre de PI) [rad]
/* -------- Electrolyte ------- */
double mu = 1e-3; // Viscosité dynamique [Pa.s]
double gravity = 9.81; // Gravity [m/s2]
double rho = 997; // Density [kg/m3]
double P0 = 1.01e5; // Pressure [Pa]
double T0 = 273.15 + 35.0; // Temperature [K]
double gammaLiquidGas[2] = {70e-3, 70e-3}; // Surface tensions [N/m]
/* ------ Electrochemistry ---- */
int Nelectronic[2] = {4, 2}; // Nombre d'électron nécessaire à la formation d'une molécule
double sigma = 2; // Electrolyte conductivity [S/m]
double tafelA0[2] = {(60e-3)/log(10), (30e-3)/log(10)}; // Tafel modeling : voltage slope [V]
double tafelJ0[2] = {1e-3, 1e-1}; // Tafel modeling : exchange current [A/m2]
/* ---------------------------- */
/* ------------ Coupling ---------- */
bool electricalCoupling = false; // Couplage faible du champ électrique
/* ------------------------------- */
/* ----------- Repeat y ---------- */
bool repeatYBubble = false; // Simule la répétition du motif selon y (recopie des bulles 1 fois)
/* ------------------------------- */
/* -------------- Mesh ------------ */
int meshN = 400; // Mesh size [-]
double BUBBLE_MESH_COEFF = 5.0; // Bubble meshing factor for bubble diameter [-]
/* -------------------------------- */
/* ------ Solver tolerance ----- */
#define TOLERANCE_LU 1e-10 // Tolérance sur la vérification de l'inverse obtenu via la méthode LU
#define TOLERANCE_POTENTIAL 1e-6 // Critère d'arrêt sur la résolution du potentiel électrolytique
#define TOLERANCE_CURRENT 1e-3 // Critère d'arrêt sur la résolution du courant
#define DERIVATIVE_VARIATION 1e-2 // Fraction de la grandeur pour l'approximation de la dérivée
#define SOLVER_MAX_ITERATION 40 // Nombre maximum d'itération (pour les cas non solvables)
/* ----------------------------- */
/*---------- Time step -------- */
double GROWING_TIME_COEF = 5.0; // Growing time step factor [-]
double BUOYANCY_TIME_COEF = 10.0; // Buoyancy time step factor [-]
double ELECTRICAL_TIME_COEF = 1.0; // Buoyancy time step factor [-]
/* ---------------------------- */
/* -- Numerical parameters -- */
/* -------------------------- */
double growingKineticFraction = 0.7; // Proportion de croissance/attente sur un site (croissance)
double detachementMinimalRadiusNorm = 0.8; // Fraction du rayon de détachement naturel à partir duquel le détachement dynamique est possible
/* -------------------------- */
#define interactionLongMaximumRange 15 // Portée maximale de l'interaction champ lointain (normalisée par le rayon moyen)
#define interactionLubMaximumRange 3.5 // Portée maximale de l'interaction lubrification (normalisée par le rayon moyen)
/*--------- Repulsive ------- */
double interactionRepulsiveRange = 2.2; // Portée du modèle de répulsion (Rref)
double interactionRepulsiveVelocity = 5e0; // Raideur du modèle de répulsion (Uref)
/* ----- Anti-collision ----- */
#define interactionCollisionRange 3 // Portée à partir de laquelle le contact est testé (normalisée par le rayon moyen)
#define interactionSafetyD 0.1 // Distance de sécurité à partir de laquelle la vitesse normale relative est absorbée
#define extraRFritz 0.02 // Rayon supplémentaire au rayon de Fritz pour éviter le blocage numérique
/* -------------------------- */
#define cloggedTimeFactor 1 // Temps à partir duquel une bulle est considérée coincée
/* -------------------------- */
/* ---------- Ouptut --------- */
/* --------------------------- */
bool outputElectrical = true;
string outputFolder = "output";
int outputMode = 0;
// Mode 0 : all frames
// Mode 1 : real time recording during the entire simulation
int outputFps = 60, outputTimeFactor = 1/10;
// Mode 2 : all frame between given interval
double outputStart = 0, outputStop = 1;
/* --------------------------- */
bool useSubProgressBar = true;
/* --------------------------- */

63
src/Simulateur/solver/Counter.h Normal file
View File

@ -0,0 +1,63 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
/* --- Compute the extraction proportions --- */
/* ------------------------------------------ */
void updateRemaining(ElectrolysisStatistics& Statistics, ListBubbles& Bubbles) {
Statistics.remaining[0] = 0; Statistics.remaining[1] = 0;
for (list<Bubble>::iterator it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) {
Statistics.remaining[0] += it->quantity[0];
Statistics.remaining[1] += it->quantity[1];
}
for (list<Bubble>::iterator it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) {
Statistics.remaining[0] += it->quantity[0];
Statistics.remaining[1] += it->quantity[1];
}
}
/* ------------------------------------------ */
/* ---------- Compute performances ---------- */
/* ------------------------------------------ */
Performances zeroPerformances() {
return {{1e9, 1e9}, {0.0, 0.0}, {0.0, 0.0}, {1e9, 1e9}, {0.0, 0.0}, 1e9};
}
Performances computePerformances(Parameters& parameters, ElectrolysisStatistics& Statistics, ListBubbles& Bubbles) {
updateRemaining(Statistics, Bubbles);
Performances performances = zeroPerformances();
for (int s = 0; s < 2; s++) {
performances.electrolysed[s] = (parameters.I*parameters.tf) / QnPrefactor[s];
performances.collected[s] = Statistics.collected[s];
performances.contamination[s] = Statistics.contamination[s];
performances.lost[s] = Statistics.lost[s];
performances.remaining[s] = Statistics.remaining[s];
}
performances.resistor = Statistics.resistor;
return performances;
}
void updateStatistics(ElectrolysisStatistics& Statistics, MatrixTension& TENSION, MatrixXd& CURRENT) {
Statistics.resistor = TENSION.resistor;
Statistics.overpotential = (TENSION.phi_anode-TENSION.phi_cathode);
Statistics.tension[0] = TENSION.min; Statistics.tension[1] = TENSION.max;
Statistics.current[0] = CURRENT.minCoeff(); Statistics.current[1] = CURRENT.maxCoeff();
}
/* ------------------------------------------ */

106
src/Simulateur/solver/Current.h Normal file
View File

@ -0,0 +1,106 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
// Calcul la carte des densités de courant sqrt(Jx2 + Jy2)
// J = sigma.grad(phi)
MatrixXd computeCurrent(Parameters& parameters, MatrixXi& MASK, MatrixTension& TENSION) {
int* N = parameters.N;
double dxInvSigma = sigma/parameters.dL[0], dyInvSigma = sigma/parameters.dL[1];
MatrixXd CURRENT(N[0], N[1]);
for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) {
// Electrolyte
if (MASK(i,j) == 0) {
double Jx = 0, Jy = 0;
// dx
if (i > 0) {
// electrolyte
if (MASK(i-1,j) == 0) Jx = dxInvSigma*(TENSION.VALUES(i,j) - TENSION.VALUES(i-1,j));
// anode
else if (MASK(i-1,j) == 1) Jx = TafelAnodeJ(TENSION.phi_anode-TENSION.VALUES(i,j));
// cathode
else if (MASK(i-1,j) == -1) Jx = TafelCathodeJ(TENSION.VALUES(i,j)-TENSION.phi_cathode);
}
// dy
int jm = (j > 0) ? j-1 : (N[1]-1);
// electrolyte
if (MASK(i,jm) == 0) Jy = dyInvSigma*(TENSION.VALUES(i,j) - TENSION.VALUES(i,jm));
// anode
else if (MASK(i,jm) == 1) Jy = TafelAnodeJ(TENSION.phi_anode-TENSION.VALUES(i,j));
// cathode
else if (MASK(i,jm) == -1) Jy = TafelCathodeJ(TENSION.VALUES(i,j)-TENSION.phi_cathode);
// norm
CURRENT(i,j) = sqrt(SQ(Jx) + SQ(Jy));
}
// Insulation
else if (MASK(i, j) == -2) {
CURRENT(i, j) = -1;
continue;
}
// Electrode (anode/cathode)
else if ((MASK(i, j) == 1) || (MASK(i, j) == -1)) {
CURRENT(i, j) = 0;
}
}
return CURRENT;
}
// Calcul la quantité de courant arrivant en chaque noeud frontière d'une électrode
// I = S.Jtafel
MatrixXd computeSourceNormalized(Parameters& parameters, MatrixXi& MASK, MatrixTension& TENSION) {
int* N = parameters.N;
double Sx = eThick*parameters.dL[1], Sy = eThick*parameters.dL[0];
MatrixXd SOURCES(N[0], N[1]);
for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) {
double src = 0;
if (MASK(i,j) == 0) {
// i-1
if (i > 0) {
// anode
if (MASK(i-1,j) == 1) src += Sx*TafelAnodeJ(TENSION.phi_anode-TENSION.VALUES(i,j));
// cathode
else if (MASK(i-1,j) == -1) src += Sx*TafelCathodeJ(TENSION.VALUES(i,j)-TENSION.phi_cathode);
}
// i+1
if (i < (N[0]-1)) {
// anode
if (MASK(i+1,j) == 1) src += Sx*TafelAnodeJ(TENSION.phi_anode-TENSION.VALUES(i,j));
// cathode
else if (MASK(i+1,j) == -1) src += Sx*TafelCathodeJ(TENSION.VALUES(i,j)-TENSION.phi_cathode);
}
// j-1
int jm = (j > 0) ? j-1 : (N[1]-1);
// anode
if (MASK(i,jm) == 1) src += Sy*TafelAnodeJ(TENSION.phi_anode-TENSION.VALUES(i,j));
// cathode
else if (MASK(i,jm) == -1) src += Sy*TafelCathodeJ(TENSION.VALUES(i,j)-TENSION.phi_cathode);
// j+1
int jp = (j < (N[1]-1)) ? j+1 : 0;
// anode
if (MASK(i,jp) == 1) src += Sy*TafelAnodeJ(TENSION.phi_anode-TENSION.VALUES(i,j));
// cathode
else if (MASK(i,jp) == -1) src += Sy*TafelCathodeJ(TENSION.VALUES(i,j)-TENSION.phi_cathode);
}
SOURCES(i,j) = src;
}
return SOURCES;
}

422
src/Simulateur/solver/Dynamics.h Normal file
View File

@ -0,0 +1,422 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
/* -------- Compute bubble evolution --------- */
/* ------------------------------------------- */
#include "utils/funinter.h"
#include "utils/fundyn.h"
#include "utils/funvelo.h"
void computeBubble(Parameters& parameters, ListBubbles& Bubbles, vector<Wall>& Walls, ElectrolysisStatistics& Statistics, double t, double dt, bool* stopFlag) {
// 0| Prepare variables
Bubbles.maxRadius = RFritzMax; Bubbles.maxVelocity = 0;
list<Bubble>::iterator it;
for (it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) {
it->interactor.count = 0; // Reset bubble-bubble interactions
}
for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) {
it->interactor.count = 0; // Reset bubble-bubble interactions
}
// 1-2| Determine interactions and build mobility matrix
std::vector<Trip> M_coefs; // Matrix coefficients for velocity system to solve
double Mij[4] { 0, 0, 0, 0 }; // Temp matrix coefficients
int i = -1;
for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) {
i++;
// 1| Test which wall and bubbles are supposed to interact significantly
checkWalls(it->interactor.walls, Walls, it->x, it->y, 0.5*interactionLongMaximumRange*it->r);
checkBubbles(it, i, Bubbles, parameters.L[1], stopFlag);
// 2| Computing mobility matrix - line i
setMobilityMatrixSphere_I(Mij);
addMobilityMatrixPlane(Mij, &(*it));
M_coefs.push_back(Trip(2*i+0, 2*i+0, Mij[0])); M_coefs.push_back(Trip(2*i+0, 2*i+1, Mij[1]));
M_coefs.push_back(Trip(2*i+1, 2*i+0, Mij[2])); M_coefs.push_back(Trip(2*i+1, 2*i+1, Mij[3]));
double sqI = SQ(it->r);
for (int jj = 0; jj < it->interactor.count; jj++) {
int j = it->interactor.interactionIndexs[jj];
if (j == -1) continue;
Bubble* jt = it->interactor.bubbles[jj];
if (jt->r == 0) continue;
double sqJ = SQ(jt->r);
int jMode = it->interactor.imagedModes[jj];
if (it->interactor.distancesSqNorm[jj] < interactionLongMaximumRangeSq) {
setMobilityMatrixSphere_J(Mij, 0.75*it->r, sqI, sqJ, it->x-jt->x, it->y-yRepeat(jt->y, parameters.L[1], jMode));
M_coefs.push_back(Trip(2*i+0, 2*j+0, Mij[0])); M_coefs.push_back(Trip(2*i+0, 2*j+1, Mij[1]));
M_coefs.push_back(Trip(2*i+1, 2*j+0, Mij[2])); M_coefs.push_back(Trip(2*i+1, 2*j+1, Mij[3]));
}
}
}
// 3| Invert mobility matrix
SpMat Minv;
if (Bubbles.countDetached > 0) {
SpMat M(2*Bubbles.countDetached, 2*Bubbles.countDetached);
M.setFromTriplets(M_coefs.begin(), M_coefs.end());
if (invertLU(Minv, M) != 0) {
// Decomposition failed
errorProgress("mobility matrix inverting error @t=" + std::to_string(t) + "s");
*stopFlag = true;
return;
}
}
// 4| Build resistance matrix and source term
std::vector<Trip> R_coefs; // Matrix coefficients for velocity system to solve
double Rij[4] { 0, 0, 0, 0 }; // Temp matrix coefficients
Eigen::VectorXd U0(2*Bubbles.countDetached); // Source term
i = -1;
for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) {
i++;
// Build source term
U0[2*i+0] = 0;
U0[2*i+1] = -bPrefactor*SQ(it->r);
// Add repulsive term
addSphereSphereRepulsive(U0, &(*it), i);
// Add inverted mobility matrix
Rij[0] = -Minv.coeff(2*i+0, 2*i+0); Rij[1] = -Minv.coeff(2*i+0, 2*i+1);
Rij[2] = -Minv.coeff(2*i+1, 2*i+0); Rij[3] = -Minv.coeff(2*i+1, 2*i+1);
// Add sphere-sphere lubrification
addSphereSphereContact_I(Rij, parameters.L[1], &(*it));
// Add sphere-plan lubrification
addSpherePlanContact_I(Rij, &(*it));
// Build resistance matrix ii
R_coefs.push_back(Trip(2*i+0, 2*i+0, Rij[0])); R_coefs.push_back(Trip(2*i+0, 2*i+1, Rij[1]));
R_coefs.push_back(Trip(2*i+1, 2*i+0, Rij[2])); R_coefs.push_back(Trip(2*i+1, 2*i+1, Rij[3]));
for (int jj = 0; jj < it->interactor.count; jj++) {
int j = it->interactor.interactionIndexs[jj];
if (j == -1) continue;
Bubble* jt = it->interactor.bubbles[jj];
if (jt->r == 0) continue;
int jMode = it->interactor.imagedModes[jj];
if (it->interactor.distancesSqNorm[jj] < interactionLongMaximumRangeSq) {
// Add inverted mobility matrix
double ratio_ji = jt->r/it->r;
Rij[0] = -ratio_ji*Minv.coeff(2*i+0, 2*j+0); Rij[1] = -ratio_ji*Minv.coeff(2*i+0, 2*j+1);
Rij[2] = -ratio_ji*Minv.coeff(2*i+1, 2*j+0); Rij[3] = -ratio_ji*Minv.coeff(2*i+1, 2*j+1);
// Add sphere-sphere lubrification
if (it->interactor.distancesSqNorm[jj] < interactionLubMaximumRangeSq) {
addSphereSphereContact_J(Rij, parameters.L[1], &(*it), jt, j, jMode);
}
// Build resistance matrix ij
R_coefs.push_back(Trip(2*i+0, 2*j+0, Rij[0])); R_coefs.push_back(Trip(2*i+0, 2*j+1, Rij[1]));
R_coefs.push_back(Trip(2*i+1, 2*j+0, Rij[2])); R_coefs.push_back(Trip(2*i+1, 2*j+1, Rij[3]));
}
}
}
// 5| Solve velocity
Eigen::VectorXd U;
if (Bubbles.countDetached > 0) {
SpMat R(2*Bubbles.countDetached, 2*Bubbles.countDetached);
R.setFromTriplets(R_coefs.begin(), R_coefs.end());
Eigen::SparseLU<SpMat> solver;
solver.compute(R);
if (solver.info() != Eigen::Success) {
// Decomposition failed
errorProgress("resistance matrix solving error @t=" + std::to_string(t) + "s");
*stopFlag = true;
return;
}
U = solver.solve(U0);
}
// 6| Update velocity of detached bubbles
i = -1;
for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) {
i++;
it->ux = U(2*i+0);
it->uy = U(2*i+1);
}
// 7| Check attached bubble detachment
double Fdyn[2];
for (it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) {
// Si le détachement dynamique est possible
if (it->r != 0) {
// Calcul des actions hydrodynamiques
Fdyn[0] = 0; Fdyn[1] = 0;
for (int kk = 0; kk < it->interactor.count; kk++) {
Bubble *kt = it->interactor.bubbles[kk];
int kMode = it->interactor.imagedModes[kk];
if (it->interactor.distancesSqNorm[kk] < interactionCollisionRangeSq) {
double deltaX = (it->x - kt->x), deltaY = (it->y - yRepeat(kt->y, parameters.L[1], kMode));
double norm = kt->ux*deltaX + kt->uy*deltaY;
if (norm > 0) {
// Si contact avec une autre bulle alors elle détache immédiatement
it->statLife = 0;
break;
}
} else if (it->interactor.distancesSqNorm[kk] < interactionLubMaximumRangeSq) {
computeSphereSphereForce(Fdyn, &(*it), kt, kMode, parameters.L[1]);
}
}
// Ajout de la pesanteur et retrait de la composante reprise par la paroi
Fdyn[1] += buoyancyPrefactor*CUBE(it->r);
double proj = Fdyn[0]*it->site->nx + Fdyn[1]*it->site->ny;
if (proj < 0) {
Fdyn[0] -= proj*it->site->nx;
Fdyn[1] -= proj*it->site->ny;
}
// Test de détachement hydrodynamique
if ((SQ(Fdyn[0]) + SQ(Fdyn[1]) > FcapilaireSq[it->site->type])) {
it->statLife = 0;
}
}
}
// 8-9| Velocity correction and integration
i = -1;
for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) {
// 8| Velocity correction
i++;
bubbleVelocityCorrection(parameters.L[1], &(*it), stopFlag);
}
for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) {
// 9| Velocity integration
it->x = it->x + it->ux*dt;
it->y = it->y + it->uy*dt;
double velocity = sqrt(SQ(it->ux) + SQ(it->uy));
if (velocity > Bubbles.maxVelocity) Bubbles.maxVelocity = velocity;
// Check for clogging
if ((it->ux == 0) && (it->uy == 0)) it->tClogg += dt;
else it->tClogg = 0;
if (it->tClogg > parameters.tClogg) {
errorProgress("clogging detected @t=" + std::to_string(t) + "s");
*stopFlag = true;
return;
}
// Update wall interactions for next step (detached bubble growing)
checkWalls(it->interactor.walls, Walls, it->x, it->y, 0.5*interactionLongMaximumRange*it->r);
}
// ?| Detached bubble growing /* !!!!!!! DISABLED (\/ see end of this file \/) !!!!!!! */
// 10| Attached bubble growing
for (it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) {
// Prise en compte du dr propre au site
Site* site = it->site;
double dr = MIN(site->growth*(dt+site->tmiss), MAX(0, RFritzExtra[site->type] - it->r));
site->tmiss = 0;
double maxDr = INF;
double itx = it->x + interactionSafetyHalfRange*dr*site->nx, ity = it->y + interactionSafetyHalfRange*dr*site->ny;
// Detection des parois ...
checkWalls(it->interactor.walls, Walls, itx, ity, 0.5*interactionLongMaximumRange*it->r);
for (int kk = 0; kk < wallMaxmimumIndicators; kk++) {
WallIndicator *indicator = &it->interactor.walls[kk];
if (indicator->distance < INF) {
double deltaR = (it->r+dr);
if (SQ(indicator->distance) < interactionSafetyHalfRangeSq*SQ(deltaR)) {
double drRemaining = interactionSafetyHalfRangeInv*indicator->distance - it->r;
if (drRemaining < maxDr) {
maxDr = drRemaining;
}
}
}
}
// ... et des bulles empechant la croissance
for (int kk = 0; kk < it->interactor.count; kk++) {
Bubble* kt = it->interactor.bubbles[kk];
int kMode = it->interactor.imagedModes[kk];
double deltaX = (kt->x - itx), deltaY = (yRepeat(kt->y, parameters.L[1], kMode) - ity), deltaR = (it->r+dr + kt->r);
double dSQ = deltaX*deltaX + deltaY*deltaY;
if (dSQ < interactionSafetyHalfRangeSq*SQ(deltaR)) {
double drRemaining = interactionSafetyHalfRangeInv*sqrt(dSQ) - it->r - kt->r;
if (drRemaining < maxDr) {
maxDr = drRemaining;
}
}
}
// Si le grossissement a été limitée, il faut le reporter
if (maxDr != INF) {
double drBis = MIN(MAX(0, maxDr), dr);
site->tmiss = (dr-drBis)/site->growth;
dr = drBis;
}
// La croissance ne peut commencer qu'à partir d'un rayon mini
if (it->r == 0) {
if ((it->r+dr) < Rmini[site->type]) {
site->tmiss += dr/site->growth;
dr = 0;
}
}
// Finalement, si grossissement il y a :
if (dr != 0) {
it->quantity[site->type] += nRPrefactor*(mathCB(it->r + dr)-mathCB(it->r));
it->r = it->r + dr;
it->x = it->x + interactionSafetyHalfRange*dr*site->nx;
it->y = it->y + interactionSafetyHalfRange*dr*site->ny;
}
}
// 11| Detached bubble exit
list<list<Bubble>::iterator> toKillDetached;
for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) {
if (it->x < 0.0) { // left exit (0)
it->statLife = 0;
// Add to stats
Statistics.collected[0] += it->quantity[0];
Statistics.contamination[0] += it->quantity[1];
Statistics.lost[1] += it->quantity[1];
} else if (it->x > parameters.L[0]) { // right exit (1)
it->statLife = 0;
// Add to stats
Statistics.collected[1] += it->quantity[1];
Statistics.contamination[1] += it->quantity[0];
Statistics.lost[0] += it->quantity[0];
} else if (it->y > parameters.L[1]) { // top recopy
it->y = it->y - parameters.L[1];
it->statLife = it->statLife - 1;
// Add to stats
if (it->statLife == 0) {
Statistics.lost[0] += it->quantity[0];
Statistics.lost[1] += it->quantity[1];
}
}
if (it->statLife == 0) {
Bubbles.countDetached--;
Bubbles.flagDetached = true;
toKillDetached.push_front(it);
continue;
}
}
if (Bubbles.flagDetached) { // If a bubble dies, reorder the list
for (list<list<Bubble>::iterator>::iterator itKill = toKillDetached.begin(); itKill != toKillDetached.end(); itKill++) Bubbles.detached.erase(*itKill);
Bubbles.flagDetached = false;
}
// 12| Attached bubble detachment
list<list<Bubble>::iterator> toKillAttached;
for (it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) {
if (it->statLife == 0) {
// Changement de statut
Bubbles.flagAttached = true;
it->statLife = bubbleLifeCount;
toKillAttached.push_front(it);
Bubbles.detached.push_back(*it);
Bubbles.countAttached--; Bubbles.countDetached++;
// Prochaine nucléation
it->site->t = t + growingKineticFractionWait*it->site->T;
it->site->occuped = false;
}
}
// 13| Reorder bubble list
if (Bubbles.flagAttached) { // If a bubble dies, reorder the list
for (list<list<Bubble>::iterator>::iterator itKill = toKillAttached.begin(); itKill != toKillAttached.end(); itKill++) Bubbles.attached.erase(*itKill);
Bubbles.flagAttached = false;
}
}
/* ------------------------------ */
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DISABLED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
// ?| Detached bubble growing
/*double dr = dt*growingFactor;
int growingBubbleCount; Bubble* growingCollisionBubbles[2]; int growingCollisionBubblesMode[2]; double growingCollisionBubblesD[2];
WallIndicator* growingCollisionWall = 0; int growingWallCount;
for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) {
growingBubbleCount = 0; growingWallCount = 0;
// Detection des bulles empechant la croissance
for (int kk = 0; kk < it->interactor.count; kk++) {
if (it->interactor.distancesSqNorm[kk] < interactionCollisionRangeSq) {
Bubble* kt = it->interactor.bubbles[kk];
int kMode = it->interactor.imagedModes[kk];
double deltaX = (kt->x - it->x), deltaY = (Y_repeat(kt->y, kMode) - it->y), deltaR = (it->r+dr + kt->r);
double dSQ = deltaX*deltaX + deltaY*deltaY;
if (dSQ < interactionSafetyHalfRangeSq*SQ(deltaR)) {
if (growingBubbleCount == 2) break;
growingCollisionBubbles[growingBubbleCount] = kt;
growingCollisionBubblesD[growingBubbleCount] = sqrt(dSQ);
growingCollisionBubblesMode[growingBubbleCount] = kMode;
growingBubbleCount++;
}
}
}
// Detection des parois empechant la croissance
for (int kk = 0; kk < wallMaxmimumIndicators; kk++) {
WallIndicator *indicator = &it->interactor.walls[kk];
if (indicator->distance < (0.5*it->r*interactionCollisionRange)) {
// Recalcul de la distance exact (à cause de l'intégration)
if (indicator->distance < interactionSafetyHalfRange*(it->r+dr)) {
if (growingWallCount == 1) break;
growingCollisionWall = indicator;
growingWallCount++;
}
}
}
if ((growingBubbleCount != 0) || (growingWallCount != 0)) {
// Simple decallage si une bulle en blocage avec une autre
if ((growingBubbleCount == 1) && (growingWallCount == 0)) {
double nx = it->x - growingCollisionBubbles[0]->x, ny = it->y - Y_repeat(growingCollisionBubbles[0]->y, growingCollisionBubblesMode[0]);
double d_inv = 1.0/growingCollisionBubblesD[0];
nx = d_inv*nx; ny = d_inv*ny;
it->x = it->x + dr*nx;
it->y = it->y + dr*ny;
// Idem avec une paroi
} else if ((growingBubbleCount == 0) && (growingWallCount == 1)) {
it->x = it->x + dr*growingCollisionWall->n.x;
it->y = it->y + dr*growingCollisionWall->n.y;
// Résolution de la position si deux bulles
} else if ((growingBubbleCount == 2) && (growingWallCount == 0)) {
continue;
// Ou une bulle et une paroi
} else if ((growingBubbleCount == 1) && (growingWallCount == 1)) {
continue;
// Pas de croissance si plus de deux bulles en blocage ou une bulle et un mur
} else {
continue;
}
}
// Croissance
it->r = it->r + dr;
double quantity = quantityFromDr(it->r, it->r + dr), tension = getTension(parameters, TENSION, it->x, it->y);
it->quantity[O2] += tension*quantity; it->quantity[H2] += (1-tension)*quantity;
if (it->r > Bubbles.maxRadius)
Bubbles.maxRadius = it->r;
}*/ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DISABLED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

163
src/Simulateur/solver/Nucleation.h Normal file
View File

@ -0,0 +1,163 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
/* ------------ Boundaries ------------ */
/* ------------------------------------ */
void updateNucleationSite(Parameters& parameters, MatrixXd& Sources, vector<Site>& Sites) {
int* N = parameters.N;
// Create a clone of sourcesValue since it will be modified
MatrixXd SourcesCopy(N[0], N[1]);
for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) {
SourcesCopy(i,j) = (double) Sources(i,j);
}
// Setup delta step from hole density and check that it's not too small.
double delta = 1.0/dSite;
int di2 = round((delta / (2.0*parameters.dL[0])) - 0.5), dj2 = round((delta / (2.0*parameters.dL[1])) - 0.5); // -0.5 car conversion longueur case
di2 = MAX(1, di2); dj2 = MAX(1, dj2);
// Update all sites
for (vector<Site>::iterator it = Sites.begin(); it != Sites.end(); ++it) {
Point pos = {it->x, it->y};
int i, j; positionToIndex(N, &i, &j, pos.x/parameters.L[0], pos.y/parameters.L[1]);
double I = 0.0; int count = 0;
// Moyennage sur un rayon de 1/2d
for (int ii = -di2; ii <= +di2; ii++) for (int jj = -dj2; jj <= +dj2; jj++) {
if ((i+ii > -1) && (i+ii < N[0]) && (j+jj > -1) && (j+jj < N[1])) {
double dii = 0.0, djj = 0.0;
if (di2 > 0) dii = ((double) ii) / di2;
if (dj2 > 0) djj = ((double) jj) / dj2;
if (mathSQ(dii) + mathSQ(djj) <= 1.0) {
if (mathSign(SourcesCopy(i+ii, j+jj)) == (it->type == 0 ? +1 : -1)) {
I += ABS(SourcesCopy(i+ii, j+jj));
SourcesCopy(i+ii, j+jj) = 0;
count++;
}
}
}
}
if (I > 0.0) {
double T = Tprefactor[it->type] / I;
it->T = T;
it->growth = RFritz[it->type]/(T*growingKineticFraction);
} else it->T = INF;
}
}
vector<Site> buildNucleationSite(Parameters& parameters, Topology& topology, MatrixXd& Sources) {
int* N = parameters.N;
mathRandInit();
// Create a clone of sourcesValue since it will be modified
MatrixXd SourcesCopy(N[0], N[1]);
for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) {
SourcesCopy(i,j) = (double) Sources(i,j);
}
// Setup delta step from hole density and check that it's not too small.
double delta = 1.0/dSite;
int di2 = round((delta / (2.0*parameters.dL[0])) - 0.5), dj2 = round((delta / (2.0*parameters.dL[1])) - 0.5); // -0.5 car conversion longueur case
if ((di2 == 0) || (dj2 == 0)) {
errorProgress("mesh refinement too small for nucleation.");
exit(1);
}
di2 = MAX(1, di2); dj2 = MAX(1, dj2);
list<Site> SiteList;
// For each poly..
for (int s = 0; s < 2; s++) for (int p = 0; p < topology.polygonCount[s]; p++) {
int n = topology.polygonIndex[s][p+1] - topology.polygonIndex[s][p];
double pos = 0.5*delta;
// For each segment...
for (int seg = 0; seg < n; seg++) {
int index = topology.polygonIndex[s][p] + seg;
int next = topology.polygonIndex[s][p] + (seg+1)%n;
Point x0 = {topology.polygonDefinitions[index].x, topology.polygonDefinitions[index].y};
Point x1 = {topology.polygonDefinitions[next].x, topology.polygonDefinitions[next].y};
// Suppression des frontière
if ( ((x0.x == 0.0) && (x1.x == 0.0)) || ((x0.x == 1.0) && (x1.x == 1.0))
|| ((x0.y == 0.0) && (x1.y == 0.0)) || ((x0.y == 1.0) && (x1.y == 1.0)) )
continue;
// DeNormalization
x0.x *= parameters.L[0]; x1.x *= parameters.L[0];
x0.y *= parameters.L[1]; x1.y *= parameters.L[1];
Vect u = {x1.x - x0.x, x1.y - x0.y};
double length = vectLength(u);
u.x = u.x / length; u.y = u.y / length;
Vect n = {u.y, -u.x};
// Parcours du segmenet
while (pos <= length) {
Point current = vectForward(x0, u, pos);
int i, j; positionToIndex(N, &i, &j, current.x/parameters.L[0], current.y/parameters.L[1]);
double I = 0.0; int count = 0;
// Moyennage sur un rayon de 1/2d
for (int ii = -di2; ii <= +di2; ii++) for (int jj = -dj2; jj <= +dj2; jj++) {
if ((i+ii > -1) && (i+ii < N[0]) && (j+jj > -1) && (j+jj < N[1])) {
double dii = 0.0, djj = 0.0;
if (di2 > 0) dii = ((double) ii) / di2;
if (dj2 > 0) djj = ((double) jj) / dj2;
if (mathSQ(dii) + mathSQ(djj) <= 1.0) {
if (mathSign(SourcesCopy(i+ii, j+jj)) == (s == 0 ? +1 : -1)) {
I += ABS(SourcesCopy(i+ii, j+jj));
SourcesCopy(i+ii, j+jj) = 0;
count++;
}
}
}
}
if (I > 0.0) {
double T = Tprefactor[s] / I;
double randomT = T*growingKineticFractionWait*mathRand();
SiteList.push_back({current.x, current.y, T, RFritz[s]/(T*growingKineticFraction), n.x, n.y, s, randomT, 0.0, 0});
}
pos += delta;
}
pos -= length;
}
}
vector<Site> Sites(SiteList.begin(), SiteList.end());
return Sites;
}
/* ------------------------------------ */
/* ----------- Nucleation ------------ */
/* ----------------------------------- */
void nucleateBubble(ListBubbles& Bubbles, vector<Site>& Sites, double t) {
if (t < Bubbles.tNextNucleation) return;
Bubbles.tNextNucleation = INF;
for (vector<Site>::iterator it = Sites.begin(); it != Sites.end(); ++it) {
if (!it->occuped) if (t >= it->t) {
// Pop a new bubble !
it->occuped = true;
it->tmiss = 0;
Bubbles.attached.push_front({it->x, it->y, 0.0, 0.0, 0.0, 1, 0, true, &(*it), {}, {0.0, 0.0}, 0.0});
Bubbles.countAttached++;
}
// May be the next one ?
if (it->t < Bubbles.tNextNucleation) {
Bubbles.tNextNucleation = it->t;
}
}
}
/* ----------------------------------- */

245
src/Simulateur/solver/Tension.h Normal file
View File

@ -0,0 +1,245 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
// Compute total current (at the anode)
double computeTotalCurrent(Parameters& parameters, MatrixXi& ID, vector<ElectrodeNode>& nodes, VectorXd& phi_n, double phi_anode) {
double I = 0;
vector<ElectrodeNode>::iterator it;
for (it = nodes.begin(); it != nodes.end(); it++) {
double phi_n_id = phi_n(it->id);
if (it->nx_anode != 0) I += it->nx_anode*eThick*parameters.dL[1]*TafelAnodeJ(phi_anode-phi_n_id);
if (it->ny_anode != 0) I += it->ny_anode*eThick*parameters.dL[0]*TafelAnodeJ(phi_anode-phi_n_id);
}
return I;
}
// Non linear part of potential equation (Tafel)
VectorXd evaluateF(Parameters& parameters, MatrixXi& ID, vector<ElectrodeNode>& nodes, VectorXd& phi_n, double phi_anode, double phi_cathode) {
double dxSigmaInv = 1.0/(sigma*parameters.dL[0]), dySigmaInv = 1.0/(sigma*parameters.dL[1]);
VectorXd F_n(phi_n.size()); F_n.setZero();
vector<ElectrodeNode>::iterator it;
for (it = nodes.begin(); it != nodes.end(); it++) {
int id = it->id;
double phi_n_id = phi_n(it->id);
if (it->nx_anode != 0) F_n(id) -= it->nx_anode*it->denomInv*dxSigmaInv*TafelAnodeJ(phi_anode-phi_n_id);
if (it->ny_anode != 0) F_n(id) -= it->ny_anode*it->denomInv*dySigmaInv*TafelAnodeJ(phi_anode-phi_n_id);
if (it->nx_cathode != 0) F_n(id) -= it->nx_cathode*it->denomInv*dxSigmaInv*TafelCathodeJ(phi_n_id-phi_cathode);
if (it->ny_cathode != 0) F_n(id) -= it->ny_cathode*it->denomInv*dySigmaInv*TafelCathodeJ(phi_n_id-phi_cathode);
}
return F_n;
}
// Jacobian of the non linear part of potential equation (Tafel)
SpMat evaluateDF(Parameters& parameters, MatrixXi& ID, vector<ElectrodeNode>& nodes, VectorXd& phi_n, double phi_anode, double phi_cathode) {
double dxSigmaInv = 1.0/(sigma*parameters.dL[0]), dySigmaInv = 1.0/(sigma*parameters.dL[1]);
SpMat dF_n(phi_n.size(), phi_n.size());
vector<Trip> coeffs;
vector<ElectrodeNode>::iterator it;
for (it = nodes.begin(); it != nodes.end(); it++) {
double phi_n_id = phi_n(it->id), dF_n_id = 0;
if (it->nx_anode != 0) dF_n_id -= it->nx_anode*it->denomInv*dxSigmaInv*dTafelAnodeJ(phi_anode-phi_n_id);
if (it->ny_anode != 0) dF_n_id -= it->ny_anode*it->denomInv*dySigmaInv*dTafelAnodeJ(phi_anode-phi_n_id);
if (it->nx_cathode != 0) dF_n_id -= it->nx_cathode*it->denomInv*dxSigmaInv*dTafelCathodeJ(phi_n_id-phi_cathode);
if (it->ny_cathode != 0) dF_n_id -= it->ny_cathode*it->denomInv*dySigmaInv*dTafelCathodeJ(phi_n_id-phi_cathode);
if (dF_n_id > 0) coeffs.push_back(Trip(it->id, it->id, dF_n_id));
}
dF_n.setFromTriplets(coeffs.begin(), coeffs.end());
return dF_n;
}
// Calcul le potentiel en chaque noeuds de la solution électrolytique
// Pour la condition : anode = +phi_anode [V], cathode = -phi_cathode [V]
VectorXd subSolveTension(Parameters& parameters, SpMat& A, MatrixXi& ID, vector<ElectrodeNode>& nodes, double phi_anode, double phi_cathode, int* failFlag) {
VectorXd phi_n(A.cols()); phi_n.setZero();
VectorXd residual, delta_phi;
int iteration = 0;
double res0 = -1;
while (1) {
// Newton-Rahpson solving
VectorXd F_n = evaluateF(parameters, ID, nodes, phi_n, phi_anode, phi_cathode);
SpMat dF_n = evaluateDF(parameters, ID, nodes, phi_n, phi_anode, phi_cathode);
residual = -(A*phi_n + F_n);
// Compute and check residual
double res = residual.norm()/phi_anode;
if (res0 == -1) res0 = res;
progressRefresh2a();
if (res < TOLERANCE_POTENTIAL) break;
// Start new iteration
if (iteration > SOLVER_MAX_ITERATION) {
*failFlag = -1;
return {};
} iteration++;
// Next step
SpMat Jacobian = A + dF_n;
*failFlag = solveLU(delta_phi, Jacobian, residual);
if (*failFlag != 0) return {};
phi_n = phi_n + delta_phi;
}
return phi_n;
}
// Calcul le potentiel en chaque noeuds de la solution électrolytique
// En respectant la condition : int(J.dS) = parameters.I
MatrixTension solveTension(Parameters& parameters, MatrixXi& MASK, bool ignoreFail) {
int* N = parameters.N;
// Parcours de cellules du mask et attribution d'id seulement dans l'électrolyte
MatrixXi ID(N[0], N[1]); ID.setZero();
int n = 0;
for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) {
if (MASK(i, j) == 0) ID(i, j) = n++;
}
double ddxInv = 1.0/mathSQ(parameters.dL[0]), ddyInv = 1.0/mathSQ(parameters.dL[1]);
// Construction de la matrice des résistance (partie linéaire) et récupération des indices à l'interface electrolyte/electrode
vector<Trip> coeffs;
vector<ElectrodeNode> nodes;
for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) if (MASK(i,j) == 0) {
int id = ID(i,j);
ElectrodeNode node = {id, 0, 0, 0, 0, 0};
bool electrolyteNeighboor[4] = {false, false, false, false};
double denom = 0;
// i-1
if (i > 0) {
if (MASK(i-1, j) == 0) {
electrolyteNeighboor[0] = true;
denom += ddxInv;
} else if (MASK(i-1, j) == 1) node.nx_anode++;
else if (MASK(i-1, j) == -1) node.nx_cathode++;
}
// i+1
if (i < (N[0]-1)) {
if (MASK(i+1, j) == 0) {
electrolyteNeighboor[1] = true;
denom += ddxInv;
} else if (MASK(i+1, j) == 1) node.nx_anode++;
else if (MASK(i+1, j) == -1) node.nx_cathode++;
}
// j-1 (with repeat condition)
int jm = (j > 0) ? j-1 : (N[1]-1);
if (MASK(i, jm) == 0) {
electrolyteNeighboor[2] = true;
denom += ddyInv;
} else if (MASK(i, jm) == 1) node.ny_anode++;
else if (MASK(i, jm) == -1) node.ny_cathode++;
// j+1 (with repeat condition)
int jp = (j < (N[1]-1)) ? j+1 : 0;
if (MASK(i, jp) == 0) {
electrolyteNeighboor[3] = true;
denom += ddyInv;
} else if (MASK(i, jp) == 1) node.ny_anode++;
else if (MASK(i, jp) == -1) node.ny_cathode++;
// i,j
coeffs.push_back(Trip(id, id, 1));
double denomInv = 1.0/denom;
if (electrolyteNeighboor[0]) coeffs.push_back(Trip(id, ID(i-1,j), -ddxInv*denomInv));
if (electrolyteNeighboor[1]) coeffs.push_back(Trip(id, ID(i+1,j), -ddxInv*denomInv));
if (electrolyteNeighboor[2]) coeffs.push_back(Trip(id, ID(i,jm), -ddyInv*denomInv));
if (electrolyteNeighboor[3]) coeffs.push_back(Trip(id, ID(i,jp), -ddyInv*denomInv));
if ((node.nx_anode != 0) || (node.ny_anode != 0) || (node.nx_cathode != 0) || (node.ny_cathode != 0)) {
// Un noeud ne peut être en contact à la fois avec l'anode et la cathode
if (((node.nx_anode != 0) || (node.ny_anode != 0)) && ((node.nx_cathode != 0) || (node.ny_cathode != 0))) {
errorProgress("electrical mesh too high (node in contact with both anode and cathode)");
exit(1);
}
// Stockage du noeud frontière pour la partie non linéaire (Tafel)
node.denomInv = denomInv;
nodes.push_back(node);
}
}
// Build linear matrix of resistors
SpMat A(n, n);
A.setFromTriplets(coeffs.begin(), coeffs.end());
// Newton-Rahpson iterative solve
int failFlag = 0;
// First voltage approximation
double phi_anode_cathode = parameters.J*parameters.L[0]/sigma + tafelA0[ANODE]*log(parameters.J/tafelJ0[ANODE]) + tafelA0[CATHODE]*log(parameters.J/tafelJ0[CATHODE]);
double d_phi_anode_cathode = DERIVATIVE_VARIATION*phi_anode_cathode;
VectorXd phi;
double I = 0;
int iteration = 0;
double res0 = -1;
while (1) {
// Compute potential and current for delta_phi
phi = subSolveTension(parameters, A, ID, nodes, phi_anode_cathode, 0.0, &failFlag);
if (failFlag != 0) break;
I = computeTotalCurrent(parameters, ID, nodes, phi, phi_anode_cathode);
// Check residual and post progress
double res = mathAbs(I-parameters.I)/parameters.I;
if (res0 == -1) res0 = res;
progressRefresh2b(mathMinMax(log(res0/res)/log(res0/TOLERANCE_CURRENT), 0.0, 1.0));
if (res < TOLERANCE_CURRENT) break;
// Start new iteration
if (iteration > SOLVER_MAX_ITERATION) {
failFlag = -2;
break;
} iteration++;
// Compute derivative
phi = subSolveTension(parameters, A, ID, nodes, phi_anode_cathode-d_phi_anode_cathode, 0.0, &failFlag);
if (failFlag != 0) break;
double dI = I - computeTotalCurrent(parameters, ID, nodes, phi, phi_anode_cathode-d_phi_anode_cathode);
// Next step (Newton)
phi_anode_cathode = phi_anode_cathode - d_phi_anode_cathode*(I-parameters.I)/(dI);
}
// Return
MatrixXd TENSION(N[0], N[1]);
if (failFlag == 0) {
// Solving success : extract result
for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) {
if (MASK(i,j) == 0) TENSION(i,j) = phi(ID(i,j));
else if (MASK(i,j) == 1) TENSION(i,j) = phi.maxCoeff();
else if (MASK(i,j) == -1) TENSION(i,j) = phi.minCoeff();
else if (MASK(i,j) == -2) TENSION(i,j) = -1.0;
else if (MASK(i,j) == -3) TENSION(i,j) = phi.maxCoeff();
else if (MASK(i,j) == -4) TENSION(i,j) = phi.minCoeff();
}
// Compute resistor (for effective surface S = e x Ly)
double J = I/(eThick*parameters.L[1]);
double resistor = (phi.maxCoeff() - phi.minCoeff()) / J;
return {TENSION, phi.minCoeff(), phi.maxCoeff(), phi_anode_cathode, 0, resistor, 0, true};
} else if (ignoreFail) {
return {TENSION, 0, 0, 0, 0, 0, 0, false};
} else {
// Solving fail : stop program
errorProgress("electrical solving error (code: " + std::to_string(failFlag) + ")");
exit(1);
}
}

44
src/Simulateur/solver/Timestep.h Normal file
View File

@ -0,0 +1,44 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
TimeSteps updateTimestep(TimeSteps& timesteps, double maxRadius, double maxVelocity, vector<Site>& Sites, bool* stopFlag) {
// Nucleation growing rate
double growingFactorMax = 0;
for (vector<Site>::iterator it = Sites.begin(); it != Sites.end(); it++)
growingFactorMax = MAX(growingFactorMax, it->growth);
double TgrowthMin = RFritzMin/growingFactorMax;
double dtMaximal1 = TgrowthMin/GROWING_TIME_COEF;
// Equilibrium buoyancy velocity
double UElec = abs(bPrefactor*mathSQ(maxRadius));
double dtElec = 2*maxRadius/(ELECTRICAL_TIME_COEF*UElec);
// Check Reynolds
double U = abs(bPrefactor*mathSQ(maxRadius));
double dtMaximal2 = 2*maxRadius/(BUOYANCY_TIME_COEF*U);
double reynolds = (rho/mu)*maxVelocity*maxRadius;
if (reynolds > 10) {
errorProgress("reynolds number too high (Re=" + std::to_string(reynolds) + ")");
*stopFlag = true;
}
// Return the smallest value
double dt = MIN(dtMaximal1, dtMaximal2);
return {dt, dtElec, timesteps.i, timesteps.iElec};
}

298
src/Simulateur/solver/Topology.h Normal file
View File

@ -0,0 +1,298 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
/* --------- Mask generation ---------- */
/* ------------------------------------ */
int checkTopologyPolySub(Topology& topology, double x, double y) {
Point X = {x, y};
for (int s = 0; s < 2; s++) for (int p = 0; p < topology.polygonCount[s]; p++) {
int n = topology.polygonIndex[s][p+1] - topology.polygonIndex[s][p];
Point *polys = new Point[n];
for (int seg = 0; seg < n; seg++) {
int index = topology.polygonIndex[s][p] + seg;
polys[seg] = {topology.polygonDefinitions[index].x, topology.polygonDefinitions[index].y};
}
if (isInside(polys, n, X)) {
return (s == 0) ? 1 : -1;
}
delete [] polys;
}
return 0;
}
int checkTopologyPoly(int* N, Topology& topology, double x, double y) {
double dY = 0.66666 / (N[1] - 1);
double r = 0;
if (r == 0) r = checkTopologyPolySub(topology, x, y + dY);
if (r == 0) r = checkTopologyPolySub(topology, x, y - dY);
return r;
}
Eigen::MatrixXi generateMask(Parameters& parameters, Topology& topology) {
int* N = parameters.N;
Eigen::MatrixXi MASK(N[0], N[1]);
// Fill mask
for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) {
double x, y; indexToPosition(N, i, j, &x, &y);
MASK(i, j) = checkTopologyPoly(N, topology, x, y); // anode(1) electrolyte(0) cathode(-1) insulation(-2) anode_trapped(-3) cathode_trapped(-4)
}
return MASK;
}
void getAttachedBubbleLocation(double& x, double& y, double& r, list<Bubble>::iterator& it) {
// Calcul du rayon actuel
r = it->r + it->site->tmiss*it->site->growth;
double dLMouillage = r*sin(M_PI - theta[it->site->type]);
// Positionnement réel d'une bulle attachée
x = it->site->x + it->site->nx*r - it->site->nx*dLMouillage;
y = it->site->y + it->site->ny*r - it->site->ny*dLMouillage;
}
void subupdateMask(Parameters& parameters, MatrixXi& MASK, list<Bubble>::iterator& it, bool attached) {
int* N = parameters.N;
double* L = parameters.L;
double itx = it->x, ity = it->y, itr = it->r;
if (attached) getAttachedBubbleLocation(itx, ity, itr, it);
int i, j; positionToIndex(N, &i, &j, itx/L[0], ity/L[1]);
double x, y;
int di = round(itr/parameters.dL[0]), dj = round(itr/parameters.dL[1]);
for (int ii = -di; ii <= di; ii++) for (int jj = -dj; jj <= dj; jj++) {
indexToPosition(N, i+ii, j+jj, &x, &y);
if ((mathSQ((L[0]*x)-itx) + mathSQ((L[1]*y)-ity)) < mathSQ(itr)) {
if (isBetween(i+ii, j+jj, N)) if (MASK(i+ii, j+jj) == 0)
MASK(i+ii, j+jj) = -2;
}
}
}
MatrixXi updateMask(Parameters& parameters, MatrixXi MASK, ListBubbles& Bubbles) {
// Replace all insulated and others specials with electrolyte
int* N = parameters.N;
for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) {
if ((MASK(i, j) == -2) || (MASK(i, j) == -3) || (MASK(i, j) == -4)) {
MASK(i, j) = 0;
}
}
// Add bubbles
list<Bubble>::iterator it;
for (it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) {
subupdateMask(parameters, MASK, it, true);
}
for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) {
subupdateMask(parameters, MASK, it, false);
}
// Remove insulated areas
// Les bulles peuvent rendre certains noeuds d'électrolyte isolés (lapluspart du temps contre une électrode)
// C'est noeuds sont alors au potentiel de l'électrode en question (J=0), mais il est nécessaire de les détecter
MatrixXi VISITED(N[0], N[1]); VISITED.setZero();
for (int i = 0; i < N[0]; i++) for (int j = 0; j < N[1]; j++) {
// Vérifier qu'un noeud d'électrolyte relie bien l'anode à la cathode
if (MASK(i,j) == 0) if (VISITED(i,j) == 0) {
bool anode = false, cathode = false;
list<Coord> toCheck = {{i,j}};
VISITED(i,j) = 1;
vector<Coord> checked;
while (toCheck.size() > 0) {
// Vérifie l'élement et propage
list<Coord>::iterator coord = toCheck.begin();
int ii = coord->i, jj = coord->j;
// voisin x-
if (ii > 0) {
if (MASK(ii-1, jj) == 0) if (VISITED(ii-1, jj) == 0) {
toCheck.push_back({ii-1, jj});
VISITED(ii-1,jj) = 1;
}
if (MASK(ii-1, jj) == 1) anode = true;
if (MASK(ii-1, jj) == -1) cathode = true;
}
// voisin x+
if (ii < (N[0]-1)) {
if (MASK(ii+1, jj) == 0) if (VISITED(ii+1, jj) == 0) {
toCheck.push_back({ii+1, jj});
VISITED(ii+1,jj) = 1;
}
if (MASK(ii+1, jj) == 1) anode = true;
if (MASK(ii+1, jj) == -1) cathode = true;
}
// voisin y-
int jjm = (jj > 0) ? jj-1 : (N[1]-1);
if (MASK(ii, jjm) == 0) if (VISITED(ii, jjm) == 0) {
toCheck.push_back({ii, jjm});
VISITED(ii,jjm) = 1;
}
if (MASK(ii, jjm) == 1) anode = true;
if (MASK(ii, jjm) == -1) cathode = true;
// voisin y+
int jjp = (jj < (N[1]-1)) ? jj+1 : 0;
if (MASK(ii, jjp) == 0) if (VISITED(ii, jjp) == 0) {
toCheck.push_back({ii, jjp});
VISITED(ii,jjp) = 1;
}
if (MASK(ii, jjp) == 1) anode = true;
if (MASK(ii, jjp) == -1) cathode = true;
// valider l'élement testé
toCheck.erase(coord);
checked.push_back({ii, jj});
}
// Si le groupe ne connecte pas l'anode à la caothde...
if (!(anode && cathode)) {
for (vector<Coord>::iterator it = checked.begin(); it != checked.end(); it++) {
if (anode && (!cathode)) MASK(it->i, it->j) = -3; // Il sera au même potentiel que l'anode
else if ((!anode) && cathode) MASK(it->i, it->j) = -4; // ... que la cathode
else MASK(it->i, it->j) = -2; // ou completement indéterminé
}
}
}
}
return MASK;
}
/* ------------------------------------ */
/* ------ Generate wall list -------- */
/* ---------------------------------- */
vector<Wall> generateWalls(Parameters& parameters, Topology& topology) {
list<Wall> WallList;
for (int s = 0; s < 2; s++) for (int p = 0; p < topology.polygonCount[s]; p++) {
int n = topology.polygonIndex[s][p+1] - topology.polygonIndex[s][p];
for (int seg = 0; seg < n; seg++) {
int index1 = topology.polygonIndex[s][p] + seg;
int index2 = topology.polygonIndex[s][p] + ((seg+1)%n);
Point P1 = {parameters.L[0]*topology.polygonDefinitions[index1].x, parameters.L[1]*topology.polygonDefinitions[index1].y};
Point P2 = {parameters.L[0]*topology.polygonDefinitions[index2].x, parameters.L[1]*topology.polygonDefinitions[index2].y};
Vect u = {P2.x-P1.x, P2.y-P1.y};
double length = vectLength(u);
u.x = u.x / length; u.y = u.y / length;
Vect n = {u.y, -u.x};
WallList.push_back({{P1.x, P1.y}, {P2.x, P2.y}, u, n, length});
WallList.push_back({{P1.x, P1.y-parameters.L[1]}, {P2.x, P2.y-parameters.L[1]}, u, n, length});
WallList.push_back({{P1.x, P1.y+parameters.L[1]}, {P2.x, P2.y+parameters.L[1]}, u, n, length});
}
}
vector<Wall> Walls(WallList.begin(), WallList.end());
return Walls;
}
/* ---------------------------------- */
/* --------- Wall distance ---------- */
/* ---------------------------------- */
WallIndicator checkWall(Wall& wall, double x, double y, double maxDistance) {
// Check if is inside
Vect P1P = {x - wall.P1.x, y - wall.P1.y};
double pos = vectProject(P1P, wall.u);
if ((pos < -maxDistance) || (pos > wall.length+maxDistance))
return {{0.0, 0.0}, {0.0, 0.0}, -1.0, nullptr};
// Test wall side
double dist; Vect n;
if (pos < 0.0) {
dist = vectLength(P1P);
if (dist > maxDistance)
return {{0.0, 0.0}, {0.0, 0.0}, -1.0, nullptr};
n = {P1P.x/dist, P1P.y/dist};
return {wall.P1, n, dist, nullptr};
} else if (pos > wall.length) {
Vect P2P = {x - wall.P2.x, y - wall.P2.y};
dist = vectLength(P2P);
if (dist > maxDistance)
return {{0.0, 0.0}, {0.0, 0.0}, -1.0, nullptr};
n = {P2P.x/dist, P2P.y/dist};
return {wall.P2, n, dist, nullptr};
} else {
dist = vectProject(P1P, wall.n);
if ((dist < 0.0) || (dist > maxDistance))
return {{0.0, 0.0}, {0.0, 0.0}, -1.0, nullptr};
n = wall.n;
return {{wall.P1.x + pos*wall.u.x, wall.P1.y + pos*wall.u.y}, n, dist, &wall};
}
}
void checkWalls(WallIndicator* indicators, vector<Wall>& Walls, double x, double y, double maxDistance) {
for (int i = 0; i < wallMaxmimumIndicators; i++) {
indicators[i].P = {-1, -1};
indicators[i].distance = INF;
}//return;
for (vector<Wall>::iterator it = Walls.begin(); it != Walls.end(); it++) {
// Test
WallIndicator indicator = checkWall(*it, x, y, maxDistance);
if (indicator.distance == -1.0)
continue;
// Une parois est détectée comme à portée...
// Si c'est un point obtenu par linéarisation (&wall = nullptr) :
// il ne doit ni déjà être dans la liste ni faire partie des extremums d'une parois déjà ajoutée
if (indicator.wall == nullptr) {
bool removeFlag = false;
for (int i = 0; i < wallMaxmimumIndicators; i++) if (indicators[i].distance != INF) {
if (indicators[i].wall == nullptr) {
if (pointEquals(indicator.P, indicators[i].P)) {
removeFlag = true;
break;
}
} else {
if (pointEquals(indicator.P, indicators[i].wall->P1) || pointEquals(indicator.P, indicators[i].wall->P2)) {
removeFlag = true;
break;
}
}
}
if (removeFlag) continue; // Rejeté !
// Sinon, il convient de vérifier que les points extremums de la parois à ajouter ne font pas l'objet d'un point obtenu par linéarisation
// Cependant, il peut y avoir plusieurs fois la parois (sur-contrainte logique)
} else {
for (int i = 0; i < wallMaxmimumIndicators; i++) if (indicators[i].distance != INF) {
if (indicators[i].wall == nullptr) {
if (pointEquals(indicator.wall->P1, indicators[i].P) || pointEquals(indicator.wall->P2, indicators[i].P)) {
// Retrait et décallage
for (int j = i; j < (wallMaxmimumIndicators-1); j++) {
indicators[j].P = indicators[j+1].P;
indicators[j].n = indicators[j+1].n;
indicators[j].distance = indicators[j+1].distance;
indicators[j].wall = indicators[j+1].wall;
}
indicators[wallMaxmimumIndicators-1].P = {-1, -1};
indicators[wallMaxmimumIndicators-1].distance = INF;
i--;
}
}
}
}
// Save the result
for (int i = 0; i < wallMaxmimumIndicators; i++) {
if (indicator.distance < indicators[i].distance) {
indicators[i].P = indicator.P;
indicators[i].distance = indicator.distance;
indicators[i].n = indicator.n;
indicators[i].wall = indicator.wall;
if (i == (wallMaxmimumIndicators-1)) maxDistance = indicator.distance;
break;
}
}
}
}
/* ---------------------------------- */

182
src/Simulateur/solver/utils/fundyn.h Normal file
View File

@ -0,0 +1,182 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
/* Dynamics functions */
/*
* Sphere-sphere long-range ------------------------------------
*/
// Mobility matrix parts for sphere-sphere longe-range interactions
void setMobilityMatrixSphere_I(double* M) {
M[0] = 1; M[1] = 0;
M[2] = 0; M[3] = 1;
}
void setMobilityMatrixSphere_J(double* M, double coeff, double sqRa, double sqRb, double rxab, double ryab) {
double r_inv, r3_inv, r5_inv, part1, part2;
r_inv = 1.0/sqrt(rxab*rxab + ryab*ryab); r3_inv = r_inv*r_inv*r_inv; r5_inv = r3_inv*r_inv*r_inv;
double Sab2 = (sqRa + sqRb) / 3.0;
part1 = (r_inv + Sab2*r3_inv); part2 = (r3_inv - 3.0*Sab2*r5_inv);
M[0] = coeff*(part1 + part2*rxab*rxab); M[1] = coeff*part2*rxab*ryab;
M[2] = coeff*part2*rxab*ryab; M[3] = coeff*(part1 + part2*ryab*ryab);
}
void addMobilityMatrixPlane(double* M, Bubble* it) {
for (int kk = 0; kk < wallMaxmimumIndicators; kk++) {
if (it->interactor.walls[kk].distance < (0.5*it->r*interactionLongMaximumRange)) {
double Ra = it->r;
double da = it->interactor.walls[kk].distance;
double daInv = 1.0/da;
double nx = it->interactor.walls[kk].n.x, ny = it->interactor.walls[kk].n.y;
double part1 = -1.125*Ra*daInv;
double part2 = -0.5625*Ra*daInv;
double part3 = (part1 - part2);
M[0] -= part2 + part3*nx*nx;
M[1] -= part3*nx*ny;
M[2] -= part3*nx*ny;
M[3] -= part2 + part3*ny*ny;
}
}
}
/*
* Sphere-sphere contact ------------------------------------
*/
// Rlub matrix parts for sphere-sphere contact interactions
void matrixRlubS(double* Rlub, bool add, double Ra, double Rb, double rxab, double ryab) {
double Reqab = Ra*Rb/(Ra+Rb);
double rabSq = (rxab*rxab + ryab*ryab);
double epsad = MAX(SMALLEST, sqrt(rabSq) - Ra - Rb)/Reqab;
double rabSqinv = 1.0/rabSq, epsadInv = 1.0/epsad;
double part1 = 0.25*epsadInv;
double part2 = 0.166666667*log(epsadInv);
double part3 = (part1 - part2)*rabSqinv;
#define Bab0 part2 + part3*rxab*rxab
#define Bab1 part3*rxab*ryab
#define Bab2 part3*rxab*ryab
#define Bab3 part2 + part3*ryab*ryab
if (add) {
Rlub[0] += Bab0;
Rlub[1] += Bab1;
Rlub[2] += Bab2;
Rlub[3] += Bab3;
} else {
Rlub[0] -= Bab0;
Rlub[1] -= Bab1;
Rlub[2] -= Bab2;
Rlub[3] -= Bab3;
}
}
// Sphere-sphere contact
// i
void addSphereSphereContact_I(double* Rij, double Ly, Bubble* it) {
for (int kk = 0; kk < it->interactor.count; kk++) {
Bubble *kt = it->interactor.bubbles[kk];
if ((it->interactor.distancesSqNorm[kk] < interactionLubMaximumRangeSq) && (kt->r != 0)) {
int kMode = it->interactor.imagedModes[kk];
matrixRlubS(Rij, false, it->r, kt->r, kt->x-it->x, yRepeat(kt->y, Ly, kMode)-it->y);
}
}
}
// j
void addSphereSphereContact_J(double* Rij, double Ly, Bubble* it, Bubble* jt, int j, int jMode) {
matrixRlubS(Rij, true, it->r, jt->r, jt->x-it->x, yRepeat(jt->y, Ly, jMode)-it->y);
}
/*
* Sphere-plan contact ------------------------------------
*/
// D matrix parts for sphere-plan contact
void matrixRlubW(double* Rlub, double Ra, double ra, double nx, double ny) {
double epsa = MAX(SMALLEST, ra - Ra)/Ra;
double epsaInv = 1.0/epsa;
double part1 = -epsaInv;
double part2 = -0.53333333333*log(epsaInv);
double part3 = (part1 - part2);
Rlub[0] += part2 + part3*nx*nx;
Rlub[1] += part3*nx*ny;
Rlub[2] += part3*nx*ny;
Rlub[3] += part2 + part3*ny*ny;
}
// Sphere-plan contact
void addSpherePlanContact_I(double* Rij, Bubble* it) {
for (int kk = 0; kk < wallMaxmimumIndicators; kk++) {
if (it->interactor.walls[kk].distance < (0.5*it->r*interactionLubMaximumRange)) {
matrixRlubW(Rij, it->r, it->interactor.walls[kk].distance, it->interactor.walls[kk].n.x, it->interactor.walls[kk].n.y);
}
}
}
/*
* Repulsive terms ------------------------------------
*/
// Sphere-sphere
void addSphereSphereRepulsive(VectorXd& U0, Bubble* it, int i) {
for (int jj = 0; jj < it->interactor.count; jj++) {
Bubble *jt = it->interactor.bubbles[jj];
// Check if k is in range and valid
if ((it->interactor.distancesSqNorm[jj] < interactionRepulsiveRangeSq) && (jt->r != 0)) {
double dxij = it->x-jt->x, dyij = it->y-jt->y;
double rSq = (SQ(dxij) + SQ(dyij));
double RrefSq = SQ(interactionRepulsiveRange);
double RrirjSq = mathSQ(it->r + jt->r);
double r = sqrt(rSq);
dxij /= r; dyij /= r;
double F = it->r*interactionRepulsiveVelocity*mathSQ((RrefSq - rSq) / (RrefSq - RrirjSq));
U0[2*i+0] -= F*dxij;
U0[2*i+1] -= F*dyij;
}
}
}
/*
* Sphere-sphere force estimation ------------------------------------
*/
void computeSphereSphereForce(double* F, Bubble* it, Bubble* jt, int jMode, double Ly) {
double Rlub[4] = {0, 0, 0, 0};
matrixRlubS(Rlub, false, it->r, jt->r, jt->x-it->x, yRepeat(jt->y, Ly, jMode)-it->y);
F[0] += Rlub[0]*jt->ux + Rlub[1]*jt->uy;
F[1] += Rlub[2]*jt->ux + Rlub[3]*jt->uy;
}

89
src/Simulateur/solver/utils/funinter.h Normal file
View File

@ -0,0 +1,89 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
/* ------ Test bubble intractability ------ */
/* ---------------------------------------- */
void testInteractibilityMode(Bubble* it, Bubble* jt, int i, int j, int modeI, int modeJ, double deltaX, double deltaY, double RavvSq, double dMax, bool* stopFlag) {
if ((mathAbs(deltaX) < dMax) && (mathAbs(deltaY) < dMax)) {
double deltaSqNorm = (deltaX*deltaX + deltaY*deltaY) / RavvSq; // Distance inter-bulle normalisée par le rayon moyen
if (deltaSqNorm < interactionLongMaximumRangeSq) {
if (it->interactor.count == interactionMaxAmountAbsolute) {*stopFlag = true; return;}
it->interactor.interactionIndexs[it->interactor.count] = j;
it->interactor.distancesSqNorm[it->interactor.count] = deltaSqNorm;
it->interactor.bubbles[it->interactor.count] = jt;
it->interactor.imagedModes[it->interactor.count] = modeJ;
it->interactor.count++;
if (jt->interactor.count == interactionMaxAmountAbsolute) {*stopFlag = true; return;}
jt->interactor.interactionIndexs[jt->interactor.count] = i;
jt->interactor.distancesSqNorm[jt->interactor.count] = deltaSqNorm;
jt->interactor.bubbles[jt->interactor.count] = it;
jt->interactor.imagedModes[jt->interactor.count] = modeI;
jt->interactor.count++;
}
}
}
void checkBubbles(list<Bubble>::iterator& it, int i, ListBubbles& Bubbles, double Ly, bool* stopFlag) {
it->statCollision = 0;
int j = i; list<Bubble>::iterator itbis = it;
// Bulles détachées
for (itbis++; itbis != Bubbles.detached.end(); itbis++) {
j++;
double Ravv = 0.5*(it->r + itbis->r), RavvSq = SQ(Ravv);
double dMax = Ravv * interactionLongMaximumRange;
/* mode 0 : pas de modifications */
double deltaX = (itbis->x - it->x), deltaY = (itbis->y - it->y);
testInteractibilityMode(&*it, &*itbis, i, j, 0, 0, deltaX, deltaY, RavvSq, dMax, stopFlag);
/* mode 1 : recopie en haut */
if (itbis->y < dMax) {
deltaX = (itbis->x - it->x), deltaY = ((itbis->y + Ly) - it->y);
testInteractibilityMode(&*it, &*itbis, i, j, 2, 1, deltaX, deltaY, RavvSq, dMax, stopFlag);
}
/* mode 2 : recopie en bas */
if (itbis->y > (Ly - dMax)) {
deltaX = (itbis->x - it->x), deltaY = ((itbis->y - Ly) - it->y);
testInteractibilityMode(&*it, &*itbis, i, j, 1, 2, deltaX, deltaY, RavvSq, dMax, stopFlag);
}
}
// Bulles attachées
for (itbis = Bubbles.attached.begin(); itbis != Bubbles.attached.end(); itbis++) {
double Ravv = 0.5*(it->r + itbis->r), RavvSq = SQ(Ravv);
double dMax = Ravv * interactionLongMaximumRange;
/* mode 0 : pas de modifications */
double deltaX = (itbis->x - it->x), deltaY = (itbis->y - it->y);
testInteractibilityMode(&*it, &*itbis, i, -1, 0, 0, deltaX, deltaY, RavvSq, dMax, stopFlag);
/* mode 1 : recopie en haut */
if (itbis->y < dMax) {
deltaX = (itbis->x - it->x), deltaY = ((itbis->y + Ly) - it->y);
testInteractibilityMode(&*it, &*itbis, i, -1, 2, 1, deltaX, deltaY, RavvSq, dMax, stopFlag);
}
/* mode 2 : recopie en bas */
if (itbis->y > (Ly - dMax)) {
deltaX = (itbis->x - it->x), deltaY = ((itbis->y - Ly) - it->y);
testInteractibilityMode(&*it, &*itbis, i, -1, 1, 2, deltaX, deltaY, RavvSq, dMax, stopFlag);
}
}
}
/* Correct y for repeating symmetry */
double yRepeat(double y, double Ly, int mode) {
if (mode == 0) return y; // mode 0
else if (mode == 1) return y+Ly; // mode 1
else return y-Ly; // mode 2
}
/* ---------------------------------------- */

86
src/Simulateur/solver/utils/funvelo.h Normal file
View File

@ -0,0 +1,86 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
/* --------- Velocity correction ---------- */
/* ---------------------------------------- */
void bubbleVelocityCorrection(double Ly, Bubble* it, bool* stopFlag);
bool bubbleVelocitySub(double Ly, Bubble* it, bool mode, WallIndicator** indicatorDone, Bubble** bubbleDone, bool* stopFlag) {
// avec les murs
for (int kk = 0; kk < wallMaxmimumIndicators; kk++) {
WallIndicator *indicator = &it->interactor.walls[kk];
if (!mode) if (indicator == *indicatorDone) continue;
if (indicator->distance < (it->r*interactionSafetyHalfRange)) {
double proj0 = it->ux*indicator->n.x + it->uy*indicator->n.y;
if (proj0 < 0) { // la bulle se déplace vers le mur ?
if (mode) {
it->ux = it->ux - proj0*indicator->n.x;
it->uy = it->uy - proj0*indicator->n.y;
*indicatorDone = indicator;
return true;
} else {
// blocage définitif ! -> arrêt
it->ux = 0; it->uy = 0;
}
}
}
}
// avec les autres spheres
for (int kk = 0; kk < it->interactor.count; kk++) {
if (it->interactor.bubbles[kk]->r != 0)
if (it->interactor.distancesSqNorm[kk] < interactionSafetyRangeSq) {
Bubble* kt = it->interactor.bubbles[kk];
int kMode = it->interactor.imagedModes[kk];
if (!mode) if (kt == *bubbleDone) continue;
if ((kt->statCollision == 0) && (it->interactor.interactionIndexs[kk] != -1)) {
// si c'est une bulle détachée non traitée, il faut lui fixer sa vitesse avant
bubbleVelocityCorrection(Ly, kt, stopFlag);
}
double P21x = (it->x - kt->x), P21y = (it->y - yRepeat(kt->y, Ly, kMode));
double proj1 = it->ux*P21x + it->uy*P21y;
if (proj1 < 0) { // la bulle considérée se déplace vers une autre ?
double proj2 = 0;
if (kt->statCollision == 2) proj2 = kt->ux*P21x + kt->uy*P21y;
if (proj1 < proj2) { // le mouvement relatif est-il bien un rapprochement ?
if (mode) {
double proj = (proj1-proj2)/(P21x*P21x + P21y*P21y);
it->ux = it->ux - proj*P21x;
it->uy = it->uy - proj*P21y;
*bubbleDone = kt;
return true;
} else {
double k = MAX(0, proj2 / proj1);
it->ux = k*it->ux;
it->uy = k*it->uy;
}
}
}
}
}
return false;
}
void bubbleVelocityCorrection(double Ly, Bubble* it, bool* stopFlag) {
if (it->statCollision != 0) return;
it->statCollision = 1;
WallIndicator* indicatorDone = nullptr; Bubble* bubbleDone = nullptr;
if (bubbleVelocitySub(Ly, it, true, &indicatorDone, &bubbleDone, stopFlag)) {
bubbleVelocitySub(Ly, it, false, &indicatorDone, &bubbleDone, stopFlag);
}
it->statCollision = 2;
}
/* ---------------------------------------- */

48
src/Simulateur/utils/loader.h Normal file
View File

@ -0,0 +1,48 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
void readConfigurationTopology(Configuration& configuration, Topology& topology, double* args) {
int pos = 0, count;
/* -------- Parameters -------- */
// System dimensions [m]
configuration.L[0] = args[pos++];
configuration.L[1] = args[pos++];
// Current density [A/m2]
configuration.J = args[pos++];
// Simulation duration [s]
configuration.tf = args[pos++];
/* --------- Topology --------- */
// Polygon count
topology.polygonCount[0] = args[pos++];
topology.polygonCount[1] = args[pos++];
// Polygon indexs
count = topology.polygonCount[0]+1;
topology.polygonIndex[0] = new int[count];
for (int i = 0; i < count; i++) {
topology.polygonIndex[0][i] = args[pos++];
}
count = topology.polygonCount[1]+1;
topology.polygonIndex[1] = new int[count];
for (int i = 0; i < count; i++) {
topology.polygonIndex[1][i] = args[pos++];
}
// Polygon defs
count = topology.polygonIndex[1][topology.polygonCount[1]];
topology.polygonDefinitions = new Point[count];
for (int i = 0; i < count; i++) {
topology.polygonDefinitions[i] = {args[pos++], args[pos++]};
}
}

174
src/Simulateur/utils/math2.h Normal file
View File

@ -0,0 +1,174 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
#define MAX(a,b) ((a)>(b) ? (a) : (b))
#define MIN(a,b) ((a)<(b) ? (a) : (b))
#define ABS(x) ((x) < 0.0 ? (-(x)) : (x))
#define SQ(x) ((x)*(x))
#define CUBE(x) ((x)*(x)*(x))
#define SMALLEST 5.0e-6
#define INF 100000.0
#define INF_NORMALIZED 100.0
/* ------ Basics ------ */
/* -------------------- */
double mathAbs(double x) {
return (x < 0.0) ? -x : x;
}
double mathSQ(double x) {
return x*x;
}
double mathCB(double x) {
return x*x*x;
}
int mathSign(double x) {
return (x == 0) ? 0 : ((x > 0) ? +1 : -1);
}
double mathMax(double a, double b) {
return (a > b) ? a : b;
}
double mathMin(double a, double b) {
return (a < b) ? a : b;
}
double mathMinMax(double x, double min, double max) {
return (x < min) ? min : ((x > max) ? max : x);
}
double mathMap(double x, double xMin, double xMax, double yMin, double yMax) {
return yMin + (x-xMin)*(yMax-yMin)/(xMax-xMin);
}
/* -------------------- */
/* ------ Random ------ */
/* -------------------- */
bool mathRandInitialized = false;
void mathRandInit() {
if (!mathRandInitialized) {
srand(0); // Constant seed for easy debug
mathRandInitialized = true;
}
}
double mathRand() {
return ((double) rand() / (RAND_MAX));
}
/* -------------------- */
/* ---- Vectors and co ---- */
/* ------------------------ */
bool pointEquals(Point& P1, Point& P2) {
return ((P1).x == (P2).x) && ((P1).y == (P2).y);
}
double vectLength(Vect& u) {
return sqrt(mathSQ(u.x) + mathSQ(u.y));
}
double vectProject(Vect& u, Vect& v) {
return u.x*v.x + u.y*v.y;
}
Point vectForward(Point& x0, Vect& u, double pos) {
return {x0.x + (u.x*pos), x0.y + (u.y*pos)};
}
/* ------------------------ */
/* ---- Index <-> position ---- */
/* ---------------------------- */
void indexToPosition(int* N, int i, int j, double *x, double *y) {
*x = ((double) i)/(N[0]-1);
*y = ((double) j)/(N[1]-1);
}
void positionToIndex(int* N, int *i, int *j, double x, double y) {
*i = round((N[0]-1)*x);
*j = round((N[1]-1)*y);
}
int indexToK(int* N, int i, int j) {
return i + (N[0]*j);
}
/* ---------------------------- */
/* --- Linear system solver --- */
/* ---------------------------- */
int solveLU(VectorXd& x, SpMat& A, VectorXd& b) {
// Decompose LU
Eigen::SparseLU<SpMat> solver;
solver.compute(A);
// Decomposition failed
if(solver.info() != Eigen::Success) return 1;
// Solve
x = solver.solve(b);
// Solving failed
if(solver.info() != Eigen::Success) return 2;
// Check tolerance
double residual = (((A*x)-b).norm() / b.norm());
if (residual > TOLERANCE_LU) return 3;
return 0;
}
int invertLU(SpMat& Ainv, SpMat& A) {
// Decompose LU
Eigen::SparseLU<SpMat> solver;
solver.compute(A);
// Decomposition failed
if(solver.info() != Eigen::Success) return 1;
// Solve
SpMat I = Eigen::MatrixXd::Identity(A.rows(), A.cols()).sparseView();
Ainv = solver.solve(I);
// Solving failed
if(solver.info() != Eigen::Success) return 2;
// Check tolerance
double residual = (A*Ainv-I).norm() / I.norm();
if (residual > TOLERANCE_LU) return 3;
return 0;
}
/* ---------------------------- */
/* ---------- Others ---------- */
/* ---------------------------- */
double bounds(double x, double min, double max) {
double range = max - min;
while (x >= max)
x -= range;
while (x < min)
x += range;
return x;
}
int isApprox(double a, double b, double tol) {
return ((a > b*(1.0-tol)) && (a < b*(1.0+tol)));
}
bool isBetween(int i, int iMin, int iMax) {
return (i >= iMin) && (i <= iMax);
}
bool isBetween(int i, int j, int* N) {
return isBetween(i, 0, N[0]-1) && isBetween(j, 0, N[1]-1);
}
/* ---------------------------- */

309
src/Simulateur/utils/output.h Normal file
View File

@ -0,0 +1,309 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
// Delete output content
namespace fs = std::filesystem;
void outputDelete(const char *folder) {
for (const auto& entry : fs::directory_iterator(folder)) {
if (entry.is_regular_file()) {
fs::remove(entry.path());
} else if (entry.is_directory()) {
outputDelete(entry.path().string().c_str());
fs::remove(entry.path());
}
}
}
// Create folder before calling loadBMP API
unsigned int dispatch_loadbmp_encode_file(const char *filename, const unsigned char *imageData, unsigned int width, unsigned int height, unsigned int components) {
fs::path filePath(filename);
string parentFolder = filePath.parent_path().string();
if (!fs::exists(parentFolder)) fs::create_directory(parentFolder);
return loadbmp_encode_file(filename, imageData, width, height, components);
}
// Text output
void outputData(const char *filename, Parameters& parameters, MatrixXd& MATRIX, double minVal, double maxVal) {
ofstream frameFile;
frameFile = ofstream(filename);
if (!frameFile.is_open()) {
cerr << "Unable to create outputfile." << endl;
exit(1);
}
frameFile << "Data:" << endl;
int width = MATRIX.rows(), height = MATRIX.cols();
double x, y;
for (int j = 0; j < (int) height; j++) for (int i = 0; i < (int) width; i++) {
int jj = height - j - 1;
indexToPosition(parameters.N, i, jj, &x, &y);
double value = MATRIX(i,j);
frameFile << parameters.L[0]*x << "," << parameters.L[1]*y << ",";
if ((value == minVal) || (value == maxVal)) frameFile << "NaN";
else frameFile << value;
frameFile << endl;
}
}
void outputData(string folder, const char *subfilename, Parameters& parameters, MatrixXd& MATRIX, double minVal, double maxVal) {
string filename = folder + string(subfilename);
outputData(filename.c_str(), parameters, MATRIX, minVal, maxVal);
}
// Bitmap output
void outputMatrix(const char *filename, MatrixXd& MATRIX, double minVal, double maxVal, double nullVal) {
int width = MATRIX.rows(), height = MATRIX.cols();
unsigned char *pixels = (unsigned char *) calloc(3*width*height, sizeof (unsigned char));
unsigned int err;
for (int i = 0; i < (int) width; i++) for (int j = 0; j < (int) height; j++) {
int jj = height - j - 1;
double v = (MATRIX(i,j) - minVal)/(maxVal-minVal);
if (MATRIX(i,j) == nullVal) {
pixels[3*(i+width*jj)+0] = 34;
pixels[3*(i+width*jj)+1] = 177;
pixels[3*(i+width*jj)+2] = 76;
} else if (v <= 0.0) {
pixels[3*(i+width*jj)+0] = 255;
pixels[3*(i+width*jj)+1] = 255;
pixels[3*(i+width*jj)+2] = 255;
} else if (v >= 1.0) {
pixels[3*(i+width*jj)+0] = 255;
pixels[3*(i+width*jj)+1] = 255;
pixels[3*(i+width*jj)+2] = 255;
} else if ((v > 0) && (v < 1.0)) {
pixels[3*(i+width*jj)+0] = mathMin(mathMax(0, mathMap(v, 0.5, 1.0, 0, 255)), 255);
pixels[3*(i+width*jj)+2] = mathMin(mathMax(0, mathMap(v, 0.0, 0.5, 255, 0)), 255);
for (int c = 0; c < 3; c++) {
pixels[3*(i+width*jj)+c] = mathMin(mathMax(0, ((int) pixels[3*(i+width*jj)+c]) + ((int) mathMin(mathMax(0, 255-mathAbs(mathMap(v, 0.25, 0.75, -255, 255))), 255))), 255);
}
}
}
err = dispatch_loadbmp_encode_file(filename, pixels, width, height, LOADBMP_RGB);
if (err) std::cout << "FAILURE: " << EXIT_FAILURE << std::endl;
free(pixels);
}
void outputMatrix(string folder, const char *subfilename, MatrixXd& MATRIX, double minVal, double maxVal, double nullVal) {
string filename = folder + string(subfilename);
outputMatrix(filename.c_str(), MATRIX, minVal, maxVal, nullVal);
}
void outputMatrix(string folder, const char *subfolder, int i, MatrixXd& MATRIX, double minVal, double maxVal, double nullVal) {
string filename = folder + string(subfolder) + string("frame_") + std::to_string(i) + string(".bmp");
outputMatrix(filename.c_str(), MATRIX, minVal, maxVal, nullVal);
}
void outputMatrixI(const char *filename, MatrixXi& MATRIX, double minVal, double maxVal) {
int width = MATRIX.rows(), height = MATRIX.cols();
unsigned char *pixels = (unsigned char *) calloc(3*width*height, sizeof (unsigned char));
unsigned int err;
for (int i = 0; i < (int) width; i++) for (int j = 0; j < (int) height; j++) {
int jj = height - j - 1;
double v = (MATRIX(i,j) - minVal)/(maxVal-minVal);
if (v == 0.0)
pixels[3*(i+width*jj)+2] = 100;
else if (v == 1.0)
pixels[3*(i+width*jj)+0] = 100;
else if ((v > 0) && (v < 1.0)) {
pixels[3*(i+width*jj)+0] = MIN(MAX(0, 255.0*v), 255);
pixels[3*(i+width*jj)+1] = MIN(MAX(0, 255.0*v), 255);
pixels[3*(i+width*jj)+2] = MIN(MAX(0, 255.0*v), 255);
} else {
pixels[3*(i+width*jj)+0] = 255;
pixels[3*(i+width*jj)+1] = 255;
pixels[3*(i+width*jj)+2] = 255;
}
}
err = dispatch_loadbmp_encode_file(filename, pixels, width, height, LOADBMP_RGB);
if (err) std::cout << "FAILURE: " << EXIT_FAILURE << std::endl;
free(pixels);
}
void outputMatrixI(string folder, const char *subfilename, MatrixXi& MATRIX, double minVal, double maxVal) {
string filename = folder + string(subfilename);
outputMatrixI(filename.c_str(), MATRIX, minVal, maxVal);
}
void outputMatrixI(string folder, const char *subfolder, int i, MatrixXi& MATRIX, double minVal, double maxVal) {
string filename = folder + string(subfolder) + string("frame_") + std::to_string(i) + string(".bmp");
outputMatrixI(filename.c_str(), MATRIX, minVal, maxVal);
}
void outputPoints(const char *filename, MatrixXd& POINTS, double normalizePoint, int RpxMax) {
int width = POINTS.rows(), height = POINTS.cols();
unsigned char *pixels = (unsigned char *) calloc(3*width*height, sizeof (unsigned char));
unsigned int err;
// Then add the points
for (int i = 0; i < (int) width; i++) for (int j = 0; j < (int) height; j++) if (POINTS(i, j) != 0.0) {
double v = MIN(MAX(-1.0, POINTS(i,j)/normalizePoint), 1.0);
int Rpx = round(RpxMax*ABS(v));
for (int ri = i-Rpx; ri <= i+Rpx; ri++) for (int rj = j-Rpx; rj <= j+Rpx; rj++) {
int rjj = height - rj - 1;
if ((ri > 0) && (ri < width-1) && (rj > 0) && (rj < height-1))
if (SQ(ri-i) + SQ(rj-j) < SQ(Rpx)) {
pixels[3*(ri+width*rjj)+0] = (v > 0.0) ? 255 : 0;
pixels[3*(ri+width*rjj)+1] = 0;
pixels[3*(ri+width*rjj)+2] = (v < 0.0) ? 255 : 0;
}
}
}
err = dispatch_loadbmp_encode_file(filename, pixels, width, height, LOADBMP_RGB);
if (err) std::cout << "FAILURE: " << EXIT_FAILURE << std::endl;
free(pixels);
}
void outputPoints(string folder, const char *subfilename, MatrixXd& POINTS, double normalizePoint, int RpxMax) {
string filename = folder + string(subfilename);
outputPoints(filename.c_str(), POINTS, normalizePoint, RpxMax);
}
void outputPoints(string folder, const char *subfolder, int i, MatrixXd& POINTS, double normalizePoint, int RpxMax) {
string filename = folder + string(subfolder) + string("frame_") + std::to_string(i) + string(".bmp");
outputPoints(filename.c_str(), POINTS, normalizePoint, RpxMax);
}
void outputPointsI(const char *filename, MatrixXi& POINTS, double normalizePoint, int RpxMax) {
int width = POINTS.rows(), height = POINTS.cols();
unsigned char *pixels = (unsigned char *) calloc(3*width*height, sizeof (unsigned char));
unsigned int err;
// Then add the points
for (int i = 0; i < (int) width; i++) for (int j = 0; j < (int) height; j++) if (POINTS(i, j) != 0.0) {
double v = MIN(MAX(-1.0, POINTS(i,j)/normalizePoint), 1.0);
int Rpx = round(RpxMax*ABS(v));
for (int ri = i-Rpx; ri <= i+Rpx; ri++) for (int rj = j-Rpx; rj <= j+Rpx; rj++) {
int rjj = height - rj - 1;
if ((ri > 0) && (ri < width-1) && (rj > 0) && (rj < height-1))
if ((ri-i)*(ri-i) + (rj-j)*(rj-j) < Rpx*Rpx) {
pixels[3*(ri+width*rjj)+0] = (v > 0.0) ? 255 : 0;
pixels[3*(ri+width*rjj)+1] = 0;
pixels[3*(ri+width*rjj)+2] = (v < 0.0) ? 255 : 0;
}
}
}
err = dispatch_loadbmp_encode_file(filename, pixels, width, height, LOADBMP_RGB);
if (err) std::cout << "FAILURE: " << EXIT_FAILURE << std::endl;
free(pixels);
}
void outputSites(const char *filename, Parameters& parameters, vector<Site>& Sites) {
MatrixXd SITES(parameters.N[0], parameters.N[1]);
for (vector<Site>::iterator it = Sites.begin(); it != Sites.end(); ++it) {
double xNorm = it->x / parameters.L[0], yNorm = it->y / parameters.L[1];
int i = MIN(MAX(0, round(parameters.N[0]*xNorm)), parameters.N[0]-1),
j = MIN(MAX(0, round(parameters.N[1]*yNorm)), parameters.N[1]-1);
SITES(i,j) = (it->type == H2 ? -1 : 1) * (1/it->T);
}
outputPoints(filename, SITES, MAX(-SITES.minCoeff(), SITES.maxCoeff()), 14);
}
void outputSites(string folder, const char *subfilename, Parameters& parameters, vector<Site>& Sites) {
string filename = folder + string(subfilename);
outputSites(filename.c_str(), parameters, Sites);
}
void outputSites(string folder, const char *subfolder, int i, Parameters& parameters, vector<Site>& Sites) {
string filename = folder + string(subfolder) + string("frame_") + std::to_string(i) + string(".bmp");
outputSites(filename.c_str(), parameters, Sites);
}
/*
* Frame saving for visualization
*/
ofstream frameFile;
bool frameFileInit = false;
double outputDt, outputNext;
void initOutputFrame(const char *filename, Parameters& parameters, TimeSteps& timesteps, double tf, Topology& topology) {
frameFile = ofstream(filename);
if (!frameFile.is_open()) {
cerr << "Unable to create outputfile." << endl;
exit(1);
}
frameFile << parameters.L[0] << " " << parameters.L[1] << " " << eThick << " " << timesteps.dtElec << " ";
int polyPos = 0;
for (int s = 0; s < 2; s++) {
frameFile << topology.polygonCount[s] << " ";
for (int i = 0; i <= topology.polygonCount[s]; i++)
frameFile << topology.polygonIndex[s][i] << " ";
int polygonPtCount = topology.polygonIndex[s][topology.polygonCount[s]];
for (int i = 0; i < polygonPtCount; i++) {
frameFile << topology.polygonDefinitions[polyPos].x << ":" << topology.polygonDefinitions[polyPos].y << " ";
polyPos++;
}
}
frameFile << "!" << endl;
if (outputMode == 1) {
outputDt = outputTimeFactor/outputFps;
outputNext = 0.0;
}
frameFileInit = true;
}
void outputFrame(string folder, const char *filename, double t, Parameters& parameters, TimeSteps& timesteps, Topology& topology, ListBubbles& Bubbles, ElectrolysisStatistics& Statistics) {
string path = folder + string(filename);
// Check init
if (!frameFileInit) initOutputFrame(path.c_str(), parameters, timesteps, parameters.tf, topology);
// Check if frame must be written
if (outputMode == 1) {
if (t < outputNext) return;
outputNext += outputDt;
} else if (outputMode == 2) {
if (t < outputStart) return;
if (t > outputStop) return;
}
// Write global data
updateRemaining(Statistics, Bubbles);
frameFile << t << " " << Statistics.resistor << ":" << Statistics.overpotential << " ";
frameFile << Statistics.tension[0] << ":" << Statistics.tension[1] << ":" << Statistics.current[0] << ":" << Statistics.current[1] << " ";
frameFile << Statistics.contamination[1] << ":" << Statistics.collected[1] << ":" << Statistics.lost[1] << ":" << Statistics.remaining[1] << " ";
// Write all bubble
list<Bubble>::iterator it;
for (it = Bubbles.attached.begin(); it != Bubbles.attached.end(); it++) {
double itx, ity, itr;
double itz = eThick/2;
getAttachedBubbleLocation(itx, ity, itr, it);
frameFile << "0:" << itx << ":" << ity << ":" << itz << ":" << itr << ":" << it->quantity[0] << ":" << it->quantity[1] << " ";
}
for (it = Bubbles.detached.begin(); it != Bubbles.detached.end(); it++) {
frameFile << "1:" << it->x << ":" << it->y << ":" << (eThick/2) << ":" << it->r << ":" << it->quantity[0] << ":" << it->quantity[1] << " ";
}
frameFile << "!" << endl;
}
void endOutputFrame(ofstream frameFile) {
frameFile.close();
}

110
src/Simulateur/utils/polygon.h Normal file
View File

@ -0,0 +1,110 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
// A C++ program to check if two given line segments intersect
// https://www.geeksforgeeks.org/how-to-check-if-a-given-point-lies-inside-a-polygon/
// Given three collinear points p, q, r, the function checks if
// point q lies on line segment 'pr'
bool onSegment(Point p, Point q, Point r) {
if (q.x <= MAX(p.x, r.x) && q.x >= MIN(p.x, r.x) &&
q.y <= MAX(p.y, r.y) && q.y >= MIN(p.y, r.y))
return true;
return false;
}
// To find orientation of ordered triplet (p, q, r).
// The function returns following values
// 0 --> p, q and r are collinear
// 1 --> Clockwise
// 2 --> Counterclockwise
int orientation(Point p, Point q, Point r) {
// See https://www.geeksforgeeks.org/orientation-3-ordered-points/
// for details of below formula.
double val = (q.y - p.y) * (r.x - q.x) -
(q.x - p.x) * (r.y - q.y);
if (val == 0) return 0; // collinear
return (val > 0)? 1 : 2; // clock or counterclock wise
}
// The main function that returns true if line segment 'p1q1'
// and 'p2q2' intersect.
bool doIntersect(Point p1, Point q1, Point p2, Point q2) {
// Find the four orientations needed for general and
// special cases
int o1 = orientation(p1, q1, p2);
int o2 = orientation(p1, q1, q2);
int o3 = orientation(p2, q2, p1);
int o4 = orientation(p2, q2, q1);
// General case
if (o1 != o2 && o3 != o4)
return true;
// Special Cases
// p1, q1 and p2 are collinear and p2 lies on segment p1q1
if (o1 == 0 && onSegment(p1, p2, q1)) return true;
// p1, q1 and q2 are collinear and q2 lies on segment p1q1
if (o2 == 0 && onSegment(p1, q2, q1)) return true;
// p2, q2 and p1 are collinear and p1 lies on segment p2q2
if (o3 == 0 && onSegment(p2, p1, q2)) return true;
// p2, q2 and q1 are collinear and q1 lies on segment p2q2
if (o4 == 0 && onSegment(p2, q1, q2)) return true;
return false; // Doesn't fall in any of the above cases
}
// Returns true if the point p lies inside the polygon[] with n vertices
bool isInside(Point polygon[], int n, Point p) {
// There must be at least 3 vertices in polygon[]
if (n < 3) return false;
// Create a point for line segment from p to infinite
Point extreme = {INF_NORMALIZED, p.y};
// Count intersections of the above line with sides of polygon
int count = 0, i = 0;
do
{
int next = (i+1)%n;
// Check if the line segment from 'p' to 'extreme' intersects
// with the line segment from 'polygon[i]' to 'polygon[next]'
if (doIntersect(polygon[i], polygon[next], p, extreme))
{
// If the point 'p' is collinear with line segment 'i-next',
// then check if it lies on segment. If it lies, return true,
// otherwise false
if (orientation(polygon[i], p, polygon[next]) == 0)
return onSegment(polygon[i], p, polygon[next]);
count++;
}
i = next;
} while (i != 0);
// Return true if count is odd, false otherwise
return count&1; // Same as (count%2 == 1)
}

126
src/Simulateur/utils/progress.h Normal file
View File

@ -0,0 +1,126 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
#define BAR_LENGTH_1 40
#define BAR_LENGTH_2 5
int currentBar1 = -1, currentBar2a = -1, currentBar2b = 0;
char temps[4] = {'|', '/', '-', '\\'};
// Init progress bar
void progressInit(JNICaller* jniCaller) {
#ifdef UPDATE_PROGRESS
UPDATE_PROGRESS(*jniCaller, 10.0);
#endif
#ifndef UPDATE_PROGRESS
string title = " Electrolysis simulation ";
int total = BAR_LENGTH_1 + 2, amount = round((total - title.length())/2);
cout << string(amount, '-') << title << string((total - amount - title.length()), '-') << endl;
cout << ".";
for (int i = 0; i < BAR_LENGTH_1; i++)
cout << " ";
cout << "." << endl << flush;
currentBar1 = 0;
currentBar2a = -1;
currentBar2b = 0;
#endif
}
// Refresh progress bar
void progressPost(int bar1, int bar2a, int bar2b, JNICaller* jniCaller) {
#ifdef UPDATE_PROGRESS
if (jniCaller != nullptr) {
float percent = round(90*t/tf) + 10.0;
if (percent != *jniCaller.percent) {
UPDATE_PROGRESS(*jniCaller, percent);
*jniCaller.percent = percent;
}
}
#endif
#ifndef UPDATE_PROGRESS
// Usage of subprogress bar
if (useSubProgressBar) {
// Refresh console only if changes
if ((currentBar1 != bar1) || (currentBar2a != bar2a) || (currentBar2b != bar2b)) {
cout << '\r';
// Bar
cout << "[";
for (int i = 0; i < BAR_LENGTH_1; i++) {
if (i < bar1) cout << "#";
else cout << " ";
}
cout << "]";
// SubBar
if (bar2b != -1) {
cout << "(";
for (int i = 0; i < BAR_LENGTH_2; i++) {
if ((i < (bar2b-1)) || (bar2b == BAR_LENGTH_2)) cout << "-";
else if (i == (bar2b-1)) cout << temps[bar2a];
else cout << " ";
}
cout << ")";
} else {
for (int i = 0; i < (2+BAR_LENGTH_2); i++) {
cout << " ";
}
}
cout << flush;
currentBar1 = bar1;
currentBar2a = bar2a;
currentBar2b = bar2b;
}
} else {
// Refresh console only if changes
if (currentBar1 != bar1) {
// Bar
if (currentBar1 == 0) cout << "[";
for (int i = 0; i < (bar1-currentBar1); i++) cout << "#";
if (bar1 == BAR_LENGTH_1) cout << "]";
cout << flush;
currentBar1 = bar1;
}
}
#endif
}
void progressRefresh1(double progress1, JNICaller* jniCaller) {
int bar1 = mathMinMax(round(progress1*BAR_LENGTH_1), 0, BAR_LENGTH_1);
progressPost(bar1, -1.0, -1.0, jniCaller);
}
void progressRefresh2b(double progress2b) {
int bar2b = round(progress2b*BAR_LENGTH_2);
progressPost(currentBar1, currentBar2a, (bar2b == 0) ? 1 : bar2b, nullptr);
}
void progressRefresh2a() {
int bar2a = (currentBar2a + 1) % 4;
if (currentBar2b == -1) progressRefresh2b(0.0);
progressPost(currentBar1, bar2a, currentBar2b, nullptr);
}
// Error
void errorProgress(const std::string &message) {
#ifndef UPDATE_PROGRESS
cout << endl << flush;
cerr << "> Simulation stopped, caused by: " << endl << message << endl << flush;
#endif
}

123
src/Simulateur/utils/typesdef.h Normal file
View File

@ -0,0 +1,123 @@
// This file is part of BUDEES, a software that allows fast simulation of bubble dynamics in electrolyzers and their interactions with electrical parameters.
//
// Copyright (C) 2024 <CNRS-ENS Rennes-Université de Rennes-CY Cergy Paris Université-Université Paris Saclay-CNAM Paris-ENS Paris-Saclay>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along with this program. If not, see https://www.gnu.org/licenses/
struct Point {
double x;
double y;
};
struct Vect {
double x;
double y;
};
struct Coord {
int i;
int j;
};
struct Topology {
int polygonCount[2];
int* polygonIndex[2];
Point* polygonDefinitions;
};
struct Configuration {
double L[2];
double J;
double tf;
};
struct Performances {
double electrolysed[2], collected[2], contamination[2], lost[2], remaining[2];
double resistor;
};
struct TimeSteps {
double dt, dtElec;
int i, iElec;
};
struct MatrixTension {
MatrixXd VALUES;
double min, max;
double phi_anode, phi_cathode;
double resistor;
double tNext;
bool sucess;
};
struct Site {
double x, y, T, growth;
double nx, ny;
int type;
double t, tmiss;
bool occuped;
};
struct Wall {
Point P1, P2;
Vect u, n;
double length;
};
#define wallMaxmimumIndicators 2
struct WallIndicator {
Point P;
Vect n;
double distance;
Wall* wall;
};
struct Bubble;
#define interactionMaxAmountAbsolute interactionLongMaximumRange*interactionLongMaximumRange
struct BubbleInteractor {
int interactionIndexs[interactionMaxAmountAbsolute];
Bubble* bubbles[interactionMaxAmountAbsolute];
double distancesSqNorm[interactionMaxAmountAbsolute];
char imagedModes[interactionMaxAmountAbsolute];
int count, countVirtual;
WallIndicator walls[wallMaxmimumIndicators];
};
struct Bubble {
double x, y, r;
double ux, uy;
int statLife, statCollision;
bool alive;
Site* site;
BubbleInteractor interactor;
double quantity[2];
double tClogg;
};
struct ListBubbles {
list<Bubble> attached;
list<Bubble> detached;
bool flagAttached, flagDetached;
int countAttached, countDetached;
double maxRadius, maxVelocity;
unsigned long tNextNucleation;
};
struct ElectrolysisStatistics {
double collected[2], lost[2], contamination[2], remaining[2];
double resistor, overpotential;
double tension[2], current[2];
};
struct ElectrodeNode {
int id;
int nx_anode, ny_anode;
int nx_cathode, ny_cathode;
double denomInv;
};