Space-dependent boundary conditions - OpenModelica

Question

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;

The figure is missing. It is often more helpful to include the Modelica code than a screenshot. — Hans Olsson, Dec 09 '22 at 16:38
Please provide enough code so others can better understand or reproduce the problem. — Community, Dec 09 '22 at 23:54

score 2 · Answer 1 · answered Dec 16 '22 at 15:39

2

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 that k_MCP will always have that value. It seems you meant Real[:,:,:] 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.

answered Dec 16 '22 at 15:39

Hans Olsson

11,123
15
38

I have added the Ailette, Glycol, Water and Octadecanol models to the problem description – Diane Jan 03 '23 at 16:31
Thank you very much for your comments. I have been able to understand my mistakes! Now I have as many variables as equations. However, there is still a problem. I update my question to include this new problem – Diane Jan 04 '23 at 09:24
When you edit your question so much good answers like this one are not relevant anymore. You should accept an answer and create a new question. That is how stackoverflow is meant to be used. We are not on a message board here. – marco Jan 04 '23 at 09:56
@marco thank you. I ask a new question on stackoverflow. Sorry it's the first time I used it! – Diane Jan 04 '23 at 10:38

Space-dependent boundary conditions - OpenModelica

1 Answers1