Diane
Reputation: 1
Space-dependent boundary conditions - OpenModelica
I want to model a storage using phase change materials (see figure for description) Basically, a hot heat transfer medium flows into one plate of the exchanger and a cold heat transfer medium flows in counterflow into another plate. The plate is located between 2 layers of encapsulated PCM. I would like to determine the outlet temperature of the hot and cold heat carriers for each instant. I made a spatial discretization (2D for the MCP part and 1D for the heat carriers).
I get more equations than variables. Below, a type of warning obtained:
[1] 11:45:13 Symbolique Erreur
Too many equations, over-determined system. The model has 3044 equation(s) and 1338 variable(s).
[2] 11:45:13 Symbolique Avertissement
[MCPBat: 1216:5-1216:43]: Equation 27 (size: 1) debit_sol = 9000.0 / (3.6 * Hot.rho_c) is not big enough to solve for enough variables.
Remaining unsolved variables are:
Already solved: Hot.rho_c
Equations used to solve those variables:
Equation 5 (size: 1): Hot.rho_c = 1000.0
This is my code:
model Octadecanol
parameter Modelica.Units.SI.Temperature T_fusion=25;
parameter Modelica.Units.SI.ThermalConductivity ks_MCP=0.301;
parameter Modelica.Units.SI.ThermalConductivity kl_MCP=0.205;
parameter Modelica.Units.SI.Density rho_MCP=850;
parameter Modelica.Units.SI.SpecificHeatCapacity cp_MCP_s=1750;
parameter Modelica.Units.SI.SpecificHeatCapacity cp_MCP_l=2150;
type Heat=Real(unit="J.kg-1");
parameter Heat h_latent = 226000;
model Ailettes
parameter Modelica.Units.SI.ThermalConductivity k_ailettes=20;
parameter Modelica.Units.SI.Volume V_ailettes=0.2;
parameter Modelica.Units.SI.Density rho_ailettes=1200;
model Glycol
parameter Modelica.Units.SI.ThermalConductivity k_c=0.6;
parameter Modelica.Units.SI.DynamicViscosity mu_c=0.001139;
parameter Modelica.Units.SI.Density rho_c=1000;
parameter Modelica.Units.SI.SpecificHeatCapacity cp_c=4186;
model Water
parameter Modelica.Units.SI.ThermalConductivity k_f=0.6;
parameter Modelica.Units.SI.DynamicViscosity mu_f=0.001139;
parameter Modelica.Units.SI.Density rho_f=1000;
parameter Modelica.Units.SI.SpecificHeatCapacity cp_f=4186;
model Batt_PCM
// Déclaration des variables
Modelica.Units.SI.Mass M_MCP;
Modelica.Units.SI.Volume V_MCP;
Modelica.Units.SI.Length L_batt_tot;
Modelica.Units.SI.Length d_h_tuyau;
Modelica.Units.SI.Area S;
Modelica.Units.SI.MassFlowRate debit_PAC;
Modelica.Units.SI.MassFlowRate debit_f;
Modelica.Units.SI.MassFlowRate debit_c;
Modelica.Units.SI.CoefficientOfHeatTransfer h_ext;
Modelica.Units.SI.CoefficientOfHeatTransfer h_f;
Modelica.Units.SI.CoefficientOfHeatTransfer h_c;
Modelica.Units.SI.CoefficientOfHeatTransfer h_cf;
Modelica.Units.SI.ThermalConductivity k_eff_MCP_ref;
Modelica.Units.SI.ThermalConductivity k_MCP_ref;
Modelica.Units.SI.SpecificHeatCapacity cp_MCP_ref;
Modelica.Units.SI.Temperature T_stock;
Modelica.Units.SI.Temperature Ts_moy_c;
Modelica.Units.SI.Temperature Ts_moy_f;
Modelica.Units.SI.DimensionlessRatio ratio_a;
// Déclaration des paramètres
parameter Integer Nombre_MCP = 2;
parameter Modelica.Units.SI.Length H_MCP = 0.5;
parameter Modelica.Units.SI.Length P_MCP = 0.5;
parameter Modelica.Units.SI.Length L_MCP = 0.4;
parameter Modelica.Units.SI.Length L_tuyau = 0.05;
parameter Modelica.Units.SI.Temperature T_init = 15+273.15;
// Déclaration des paramètres de discrétisation
parameter Real Delta_x = 0.05;
parameter Real Delta_y = 0.05;
parameter Integer x_tot = integer(L_MCP/Delta_x);
parameter Integer y_tot = integer(H_MCP/Delta_y);
Integer x;
Integer y;
Integer k;
// Déclaration des constantes
final constant Real pi=2*Modelica.Math.asin(1.0); // 3.14159265358979;
// Modèles intégrant les propriétés thermophysiques des fluides caloporteurs, du MCP et des ailettes
Water Cold;
Glycol Hot;
Octadecanol PCM;
Ailettes Fin;
// Modèle intégrant le comportement de la façade solaire : donne la température de sortie de la façade solaire -> entrée du fluide caloporteur chaud
// Panneaux_Batisol Panels; // Penser à changer T_ext par Panels.T_ext
parameter Modelica.Units.SI.Temperature T_ext = 15+273.15;
parameter Modelica.Units.SI.MassFlowRate debit_sol;
parameter Modelica.Units.SI.Temperature T_out_PAC = 10+273.15;
parameter Modelica.Units.SI.Temperature T_out_sol = 100+273.15;
Real[:,:,:] k_MCP = fill(0,x_tot,y_tot,Nombre_MCP+1);
Real[:,:,:] k_eff_MCP = fill(0,x_tot,y_tot,Nombre_MCP+1);
Real[:,:,:] cp_MCP = fill(0,x_tot,y_tot,Nombre_MCP+1);
Real[:,:] Tc_y = fill(1,y_tot,Nombre_MCP) * T_init;
Real[:,:] Tf_y = fill(1,y_tot,Nombre_MCP) * T_init;
Real[:,:,:] T_xy = fill(1,x_tot,y_tot,Nombre_MCP+1)*PCM.T_fusion;
Real[:,:,:] fl_xy = fill(0,x_tot,y_tot,Nombre_MCP+1);
Real[:] Poly = fill(0,3); // Vecteur pour mettre les coefficients du polynôme de degré 2
Real[:,:] racines = fill(0,2,2); // Matrice pour mettre les racines du polynome de degré 2 (colonne 1 : partie réelle / colonne 2 : partie imaginaire)
// Coefficients calculés pour les équations du modèle de la batterie MCP
Real A, A_prim, A_scnd, B, B_prim, C, C_x_prim, C_x_scnd, C_y_prim, C_y_scnd, C_xy, C_xy_prim, C_xy_scnd, C_xy_trio, CC_x, CC_y, CC_xy, CC_xy_prim, CC_xy_scnd, CC_xy_trio, D, D_x_prim, D_x_scnd, D_y_prim, D_y_scnd, D_xy, D_xy_prim, D_xy_scnd, D_xy_trio, DD_x, DD_y, E, F, F_scnd, G, G_scnd, H, H_scnd, I, I_scnd, J_prim, J_scnd, K_prim, K_scnd, L, L_scnd, M, M_scnd, N, N_scnd, O, O_scnd;
// Real K, J, H_prim, I_prim;
initial equation
Ts_moy_c = T_init;
Ts_moy_f = T_init;
T_stock = T_init;
equation
debit_sol = 60 * 150 / (3.6*Hot.rho_c);
M_MCP = 50;
V_MCP = M_MCP / PCM.rho_MCP;
L_batt_tot = (L_MCP+2)*Nombre_MCP + L_MCP;
ratio_a = V_MCP / (V_MCP+Fin.V_ailettes);
S = P_MCP*L_tuyau;
d_h_tuyau = (2*S^2)/(P_MCP+L_tuyau);
debit_PAC = 2.5;
debit_f = debit_PAC / Nombre_MCP;
//debit_c = Panels.debit_sol / Nombre_MCP;
debit_c = debit_sol / Nombre_MCP;
h_ext = 500;
h_f = 0.664*debit_f^0.5*Cold.k_f^(2/3)*Cold.cp_f^(1/3) / (Cold.mu_f^(1/6)*d_h_tuyau);
h_c = 0.664*debit_c^0.5*Hot.k_c^(2/3)*Hot.cp_c^(1/3) / (Hot.mu_c^(1/6)*d_h_tuyau);
h_cf = 1/(1/h_c+1/h_f);
// INITIALISATION
k_MCP_ref = PCM.ks_MCP + fl_xy[1,1,1]*(PCM.kl_MCP-PCM.ks_MCP);
k_eff_MCP_ref = ratio_a*k_MCP_ref + (1-ratio_a)*Fin.k_ailettes;
cp_MCP_ref = PCM.cp_MCP_s + fl_xy[1,1,1]*(PCM.cp_MCP_l-PCM.cp_MCP_s);
// MODELISATION
// CONDITIONS AUX LIMITES (y = 1)
for y in 1:1 loop
// Initialisation : on commence par la partie de la batterie gauche (MCP + fluides chaud et froid), avant de déterminer pour les k batteries (Nombre_MCP)
k = 1;
// PARTIE MCP
// Conditions aux limites (x = 1 et y = 1)
for x in 1:1 loop
A = (ratio_a * PCM.ks_MCP + (1 - ratio_a) * Fin.k_ailettes) / (2 * Delta_y);
B = ratio_a * (PCM.kl_MCP - PCM.ks_MCP) / (2 * Delta_y);
A_prim = T_xy[x,y+2,k]-4*T_xy[x,y+1,k];
C_xy = ((-fl_xy[x,y+2,k]+4*fl_xy[x,y+1,k])/(2*Delta_y)*(-T_xy[x,y+2,k]+4*T_xy[x,y+1,k])/(2*Delta_y) + (-fl_xy[x+2,y,k]+4*fl_xy[x+1,y,k])/(2*Delta_x)*(-T_xy[x+2,y,k]+4*T_xy[x+1,y,k])/(2*Delta_x))*ratio_a*(PCM.ks_MCP-PCM.kl_MCP) - ((-5*T_xy[x,y+1,k]+4*T_xy[x,y+2,k]-T_xy[x,y+3,k])/Delta_y^2 + (-5*T_xy[x+1,y,k]+4*T_xy[x+2,y,k]-T_xy[x+3,y,k])/Delta_x^2)*(ratio_a*PCM.ks_MCP+(1-ratio_a)*Fin.k_ailettes);
CC_xy = ratio_a*(PCM.kl_MCP-PCM.ks_MCP)*((-5*T_xy[x,y+1,k]+4*T_xy[x,y+2,k]-T_xy[x,y+3,k])/Delta_y^2 + (-5*T_xy[x+1,y,k]+4*T_xy[x+2,y,k]-T_xy[x+3,y,k])/Delta_x^2) + ratio_a*(PCM.kl_MCP-PCM.ks_MCP)*((-3*T_xy[x+2,y,k]+12*T_xy[x+1,y,k])/(4*Delta_x^2) + (-3*T_xy[x,y+2,k]+12*T_xy[x,y+1,k])/(4*Delta_y^2));
D_xy = -(3*ratio_a*(PCM.kl_MCP-PCM.ks_MCP) * ((-fl_xy[x,y+2,k]+4*fl_xy[x,y+1,k])/(4*Delta_y^2) + (-fl_xy[x+2,y,k]+4*fl_xy[x+1,y,k])/(4*Delta_x^2)) + 2*(ratio_a*PCM.ks_MCP+(1-ratio_a)*Fin.k_ailettes)*(1/Delta_x^2+1/Delta_y^2));
fl_xy[x,y,k] = (PCM.rho_MCP*PCM.h_latent*der(fl_xy[x,y,k]) + C_xy + D_xy*T_xy[x,y,k]) / (CC_xy+T_xy[x,y,k]*ratio_a*(PCM.kl_MCP-PCM.ks_MCP)*(17/(4*Delta_y^2)+17/(4*Delta_x^2)));
T_xy[x,y,k] = (h_ext*T_ext + A*A_prim + B*A_prim*fl_xy[x,y,k]) / (h_ext - 3*A - 3*B*fl_xy[x,y,k]);
end for;
// Conditions aux limites (x = 2:x_tot-1 et y = 1)
for x in 2:x_tot-1 loop
C_y_prim = ((-fl_xy[x,y+2,k]+4*fl_xy[x,y+1,k])/(2*Delta_y)*(-T_xy[x,y+2,k]+4*T_xy[x,y+1,k])/(2*Delta_y) + (fl_xy[x+1,y,k]-fl_xy[x-1,y,k])/(2*Delta_x)*(T_xy[x+1,y,k]-T_xy[x-1,y,k])/(2*Delta_x))*ratio_a*(PCM.ks_MCP-PCM.kl_MCP) - ((-5*T_xy[x,y+1,k]+4*T_xy[x,y+2,k]-T_xy[x,y+3,k])/Delta_y^2 + (T_xy[x+1,y,k]+T_xy[x-1,y,k])/Delta_x^2)*(ratio_a*PCM.ks_MCP+(1-ratio_a)*Fin.k_ailettes);
CC_y = ratio_a*(PCM.kl_MCP-PCM.ks_MCP)*((-5*T_xy[x,y+1,k]+4*T_xy[x,y+2,k]-T_xy[x,y+3,k])/Delta_y^2 + (T_xy[x+1,y,k]+T_xy[x-1,y,k])/Delta_x^2) + ratio_a*(PCM.kl_MCP-PCM.ks_MCP)*(3*T_xy[x,y+2,k]-12*T_xy[x,y+1,k])/(4*Delta_y^2);
D_y_prim = -(3*ratio_a*(PCM.kl_MCP-PCM.ks_MCP) * (fl_xy[x,y+2,k]-4*fl_xy[x,y+1,k])/(4*Delta_y^2) + 2*(ratio_a*PCM.ks_MCP+(1-ratio_a)*Fin.k_ailettes)*(1/Delta_y^2-1/Delta_x^2));
fl_xy[x,y,k] = (PCM.rho_MCP*PCM.h_latent*der(fl_xy[x,y,k]) + C_y_prim + D_y_prim*T_xy[x,y,k]) / (CC_y+T_xy[x,y,k]*ratio_a*(PCM.kl_MCP-PCM.ks_MCP)*(17/(4*Delta_y^2)-2/Delta_x^2));
E = der(fl_xy[x,y,k])*(PCM.cp_MCP_s-PCM.cp_MCP_l);
Poly[1] = E*ratio_a*(PCM.kl_MCP-PCM.ks_MCP)*(17/(4*Delta_y^2)-2/Delta_x^2);
Poly[2] = E*CC_y - D_y_prim*der(T_xy[x,y,k])*(PCM.cp_MCP_l-PCM.cp_MCP_s) - E*PCM.T_fusion*ratio_a*(PCM.kl_MCP-PCM.ks_MCP)*(17/(4*Delta_y^2)-2/Delta_x^2) - PCM.cp_MCP_s*der(T_xy[x,y,k])*ratio_a*(PCM.kl_MCP-PCM.ks_MCP)*(17/(4*Delta_y^2)-2/Delta_x^2);
Poly[3] = - (E*PCM.T_fusion*CC_y - PCM.cp_MCP_s*der(T_xy[x,y,k])*CC_y - C_y_prim*der(T_xy[x,y,k])*(PCM.cp_MCP_l-PCM.cp_MCP_s) - PCM.rho_MCP*PCM.h_latent*der(fl_xy[x,y,k])*der(T_xy[x,y,k])*(PCM.cp_MCP_l-PCM.cp_MCP_s));
racines = Modelica.Math.Polynomials.roots({Poly[1],Poly[2],Poly[3]});
T_xy[x,y,k] = racines[1,1];
end for;
// Conditions aux limites (x = x_tot et y = 1)
for x in x_tot:x_tot loop
A = (ratio_a * PCM.ks_MCP + (1 - ratio_a) * Fin.k_ailettes) / (2 * Delta_y);
B = ratio_a * (PCM.ks_MCP - PCM.kl_MCP) / (2 * Delta_y);
G = T_xy[x,y+2,k]-4*T_xy[x,y+1,k];
C_xy_trio = ((-fl_xy[x,y+2,k]+4*fl_xy[x,y+1,k])/(2*Delta_y)*(-T_xy[x,y+2,k]+4*T_xy[x,y+1,k])/(2*Delta_y) + (fl_xy[x-2,y,k]-4*fl_xy[x-1,y,k])/(2*Delta_x)*(T_xy[x-2,y,k]-4*T_xy[x-1,y,k])/(2*Delta_x))*ratio_a*(PCM.ks_MCP-PCM.kl_MCP) - ((-5*T_xy[x,y+1,k]+4*T_xy[x,y+2,k]-T_xy[x,y+3,k])/Delta_y^2 + (-5*T_xy[x-1,y,k]+4*T_xy[x-2,y,k]-T_xy[x-3,y,k])/Delta_x^2)*(ratio_a*PCM.ks_MCP+(1-ratio_a)*Fin.k_ailettes);
CC_xy_trio = ratio_a*(PCM.kl_MCP-PCM.ks_MCP)*((-5*T_xy[x,y+1,k]+4*T_xy[x,y+2,k]-T_xy[x,y+3,k])/Delta_y^2 + (-5*T_xy[x-1,y,k]+4*T_xy[x-2,y,k]-T_xy[x-3,y,k])/Delta_x^2) + ratio_a*(PCM.kl_MCP-PCM.ks_MCP)*((3*T_xy[x-2,y,k]-12*T_xy[x-1,y,k])/(4*Delta_x^2) + (-3*T_xy[x,y+2,k]+12*T_xy[x,y+1,k])/(4*Delta_y^2));
D_xy_trio = -(3*ratio_a*(PCM.kl_MCP-PCM.ks_MCP) * ((-fl_xy[x,y+2,k]+4*fl_xy[x,y+1,k])/(4*Delta_y^2) + (fl_xy[x-2,y,k]-4*fl_xy[x-1,y,k])/(4*Delta_x^2)) + 2*(ratio_a*PCM.ks_MCP+(1-ratio_a)*Fin.k_ailettes)*(1/Delta_x^2+1/Delta_y^2));
fl_xy[x,y,k] = (PCM.rho_MCP*PCM.h_latent*der(fl_xy[x,y,k]) + C_xy_trio + D_xy_trio*T_xy[x,y,k]) / (CC_xy_trio+T_xy[x,y,k]*ratio_a*(PCM.kl_MCP-PCM.ks_MCP)*(17/(4*Delta_y^2)+17/(4*Delta_x^2)));
T_xy[x,y,k] = (Tc_y[y,k]*h_c + A*G + B*G*fl_xy[x,y,k]) / (h_c - 3*A - 3*B*fl_xy[x,y,k]);
end for;
// PARTIE FLUIDES CHAUD ET FROID
H = h_cf*Delta_y + h_c*Delta_y - (3*debit_c*Hot.cp_c)/(2*Delta_y);
I = -Hot.rho_c*Delta_y*L_tuyau*Hot.cp_c*der(Tc_y[y,k]) - debit_c*Hot.cp_c*(4*Tc_y[y+1,k]-Tc_y[y+2,k])/(2*Delta_y);
Tc_y[y,k] = (I+h_cf*Tf_y[y,k]+h_c*Delta_y*T_xy[x_tot,y,k]) / H;
// J = h_f*Delta_y-h_cf*Delta_y+(3*debit_f*Cold.cp_f)/(2*Delta_y);
// K = Cold.rho_f*Delta_y*L_tuyau*Cold.cp_f*der(Tf_y[y,k]) + debit_f*Cold.cp_f*(4*Tf_y[y+1,k]-Tf_y[y+2,k])/(2*Delta_y);
// Tf_y[y,k] = (K+h_f*Delta_y*T_xy[1,y,k+1]-h_cf*Delta_y*Tc_y[y,k]) / J;
Tf_y[y,k] = T_out_PAC;
end for;
end Batt_PCM;
Upvotes: 0
Views: 112
Answers (1)
Hans Olsson
Reputation: 12507
It's still not possible to fully understand the cause since the models for Water, Glycol, Octadecanol, and Ailettes are missing.
However, even without that it seems you are not using the equation-section in the right way.
There are two parts to that:
- Equation-sections aren't algorithms Difference between equation and algorithm section
- Bindings for variables are also equations, so
Real[:,:,:] k_MCP = fill(0,x_tot,y_tot,Nombre_MCP+1);means thatk_MCPwill always have that value. It seems you meantReal[:,:,:] k_MCP (start= fill(0,x_tot,y_tot,Nombre_MCP+1));at least for some of the variables.
The second item is causing the "too many equations" issue. In general, try to construct smaller models first and test them separately.
Upvotes: 2
Related Questions
- Does this OpenModelica model make sense in a physical point of view?
- Discretizing time derivative terms of PDE in the Modelica model
- Discretizing PDE in space for use with modelica
- OpenModelica solving a PDE Initialization error
- Instantiation of equation
- Singular inconsistent scalar system when connecting a pressure boundary
- Boundary assignment in openFOAM
- 1D Heat Diffusion PDE implementation in Modelica(Dymola)
- Boundary conditions in differential equations in Modelica
- Initial Conditions in OpenModelica