Почему я встретил деление на ноль ошибок в этой модели под платформой Modelica - PullRequest
Новая
       13

Почему я встретил деление на ноль ошибок в этой модели под платформой Modelica

0 голосов
/ 20 марта 2020

Я делаю модель двигателя внутреннего сгорания. Ниже приведено сообщение об ошибке, с которым я столкнулся.

dymosim запущен ... загрузка "dsin.txt" (входной файл dymosim) Во время была обнаружена следующая ошибка: 0 Ошибка: единственная несовместимая скалярная система для der (u) = (- (если Period == 61, то p-101300.0 иначе (если Period == 23 или Period == 41, то der (m) u- (der (Q) -p der) (V) + h_in der (m_in) + h_out der (m_out)) else (если Period == 12, то p - (c* (thetta-thetta_1) ^ 2 + d) else p- () / ((если Period == 61, то 0.0 else (если Period == 23 или Period == 41, то m else 0.0))) = -101300/0 Ошибка: не удалось запустить модель.

В сообщении об ошибке я мог видеть деление на ноль, но я не знаю почему, потому что я определенно установил, что начальные уравнения не делятся на нулевое деление во время = 0.

Пожалуйста, помогите мне! 

model enginemodel_ContinuousEditEq

  import SI = Modelica.SIunits;
  import nonSI = Modelica.SIunits.Conversions.NonSIunits;
  import Modelica.SIunits.Conversions.*;
  final constant Real pi=2*Modelica.Math.asin(1.0); // 3.14159265358979;

  // Medium : Air(N2 76.8%, O2 23.2%) + Fuel(C8H18 Octane 100%) => AFR 14.7 //

  replaceable package Medium_combustion = Secondyear_start.CombustionMix;
  replaceable package Medium_emission = Secondyear_start.CombustionEmission;

  CombustionMix.ThermodynamicState Mixture_st "Burned Mixture의 state";
  CombustionMix.ThermodynamicState Mixture_st_in "Burned Mixture의 inlet state";
  CombustionEmission.ThermodynamicState Emission_st "Emission Gas의 state";

  SI.Pressure p "실린더 압력 when fired";
  SI.Pressure p2 "Cycle p2";
  SI.Pressure p_c "실린더 압력 when motored";
  SI.Pressure p_wos "실린더 압력 in Woschni Equation";
  SI.Pressure p_c_wos "실린더 압력 when motored in Woschni Equation";
  SI.Volume V "실린더 체적";
  SI.Volume V_c "TDC 체적";
  SI.Volume V_s "Swept 체적";
  SI.Mass m "Mixture gas 질량";
  SI.Mass m_in "입력 질량";
  SI.Mass m_out "출력 질량";
  SI.Mass m_air "공기 질량";
  SI.Mass m_fuel "연료 질량";
  SI.SpecificInternalEnergy u "Cylinder 내부에너지";
  SI.SpecificEnthalpy h_in "Inlet 엔탈피";
  SI.SpecificEnthalpy h_out "Outlet 엔탈피";
  SI.Heat Q "Net Heat";
  SI.Heat Q_x "Combustion of the charge에 의한 발열량";
  SI.Heat Q_w "Wall로 전달된 열량";
  SI.Heat Q_t "연소 중 연료에 의한 총 발열량";
  SI.Power Q_w_dot "der(Q_w)";
  SI.Velocity w "Woschni 식 변수 중 1개";
  SI.Velocity C_m "Mean Piston Speed";
  nonSI.Angle_deg thetta "Crank Angle";
  SI.Angle thetta_rad "rad 단위 각도";
  nonSI.AngularVelocity_rpm omega "Engine RPM";
  SI.Length S_p "Swept Length of piston";
  SI.Temperature T "Gas Temp";
  Real Alpha;
  SI.CoefficientOfHeatTransfer Alpha_w;
  Real x_b;
  Real R_gas(final unit="J/(kg.K)") "Molar gas constant";
  Integer Period "Cycle 단계 flag";
  SI.Area A_w "Contact Area of the flame";
  SI.Area A_w_max "Area of the Wall";
  SI.HeatFlux WallHeatFlux "WallHeatFlux through Cylinder";

  SI.Temperature T_w "Wall Temp";
  SI.Temperature T_co "냉각수 입구 온도";
  SI.Temperature T2 "Cycle T2";
  SI.TemperatureDifference T_LMTD "대수평균온도차";
  SI.TemperatureDifference del_T1 "온도차1 for LMTD";
  SI.TemperatureDifference del_T2 "온도차2 for LMTD";
  parameter SI.TemperatureDifference beta = 0.7;
  parameter Real zeta = 15 "Penalty factor in the Robust LMTD method";
  parameter SI.TemperatureDifference tol = 0.0001;
  Real U(final unit="W/(m2.K)") "총괄 열전달 계수";
  SI.CoefficientOfHeatTransfer h_c "냉각수 대류 열전달 계수";
  SI.ReynoldsNumber Re_Dh "냉각채널 레이놀즈 수";
  SI.NusseltNumber Nu_Dh "냉각채널 누셀 수";
  SI.PrandtlNumber Pr "냉각채널 프란틀 수";
  Real f_c "friction factor by Petukhov";
  SI.Length D_h "냉각채널 수력직경";
  SI.Pressure del_p_c "냉각채널 압력 강하량";
  SI.KinematicViscosity vis_c "냉각수 동점성계수";
  SI.Velocity v_c "냉각수 유속";
  SI.VolumeFlowRate Q_c "냉각수의 volume flow rate";
  Real lpm_c "냉각수의 lpm";
  Real c "isentropic pressure parameter 1";
  Real d "isentropic pressure parameter 2";
  Real e "isentropic Temperature parameter 1";
  Real f "isentropic Temperature parameter 2";
  Integer i "iteration #";

  parameter SI.Temperature T_ci=273+25 "냉각수 입구 온도";
  parameter SI.SpecificHeatCapacityAtConstantPressure Cpc=4193 "specific heat of Coolant";
  parameter Real epsilon=8.6 "압축비";
  parameter Real LHV_f(unit="J/kg") = 44.7*(10^6) "옥탄 LHV";
  parameter Real AFR = 15;
  parameter Integer b=1 "Shape Factor";
  parameter Real x_star=0.999 "Burned fraction at the end of combustion";
  parameter nonSI.Angle_deg del_thetta = 30 "Combustion Duration";
  parameter SI.Length L=144*(10^(-3)) "커넥션 로드 길이";
  parameter SI.Length B=82.7*(10^(-3)) "실린더 보어";
  parameter SI.Length R=45*(10^(-3)) "크랭크샤프트 반지름";
  parameter SI.MassFlowRate m_dot_in=60000e-6 "혼합 공기 유량, 600mg/180deg";
  parameter nonSI.Angle_deg thetta_0=340 "Ignition Angle";
  parameter nonSI.Angle_deg thetta_1=180 "Compression Angle";
  parameter nonSI.Angle_deg thetta_3=370 "Expansion Angle";
  parameter nonSI.Angle_deg thetta_4=540 "Exhaust Angle";
  parameter SI.ThermalConductivity k_w=160 "엔진 벽 열전도계수(알루미늄 주물)";
  parameter SI.ThermalConductivity k_c=0.670 "냉각수 열전도계수";
  parameter SI.DynamicViscosity u_c=0.378*10^(-3) "Dynamic Viscosity of Coolant";
  parameter SI.Density rho_c=975 "Density of Coolant";
  parameter SI.Length L_cooling=0.2 "냉각채널 길이";
  parameter SI.Length D_c=10*0.001 "냉각채널 가로길이";
  parameter SI.MassFlowRate m_dot_c=1/4 "냉각수 유량";
  parameter SI.Length t_w=15*0.001 "Thickness of Engine wall";
  parameter Real k=1.35 "비열비";

  Modelica.Blocks.Sources.RealExpression RPM_Engine(y=3000) annotation (Placement(transformation(extent={{-14,-10},{6,10}})));

equation 
  Q_c=m_dot_c/rho_c;
  lpm_c=60000*Q_c;

  //Gas Properties//
  Mixture_st = Medium_combustion.setState_pTX(p,T,{1316/1830, 400/1830, 114/1830});
  Mixture_st_in = Medium_combustion.setState_pTX(101.3e3, 273+60,{1316/1830, 400/1830, 114/1830});
  Emission_st= Medium_emission.setState_pTX(p, T, {352/1830, 162/1830, 1316/1830});

  u = Medium_combustion.specificInternalEnergy(Mixture_st);
  R_gas = Medium_combustion.gasConstant(Mixture_st);
  h_in = Medium_combustion.specificEnthalpy(Mixture_st_in);
  h_out = Medium_emission.specificEnthalpy(Emission_st);

  RPM_Engine.y=omega;

  der(thetta)=(180/pi)*(from_rpm(omega)) "각속도와 시간 연결";

  when thetta >= 720 then
    reinit(thetta, 0);
  end when;

  thetta_rad=from_deg(thetta) "deg를 rad으로 변경";

  //크랭크 각에 따른 행정 상태 Flag(Period)//
  when 0<=thetta and thetta<180 then
    Period=61; //흡입//
  elsewhen 180<=thetta and thetta<340 then
    Period=12; //압축//
  elsewhen 340<=thetta and thetta<370 then
    Period=23; //폭발//
  elsewhen 370<=thetta and thetta<540 then
    Period=34; //팽창//
  elsewhen 540<=thetta and thetta<720 then
    Period=41; //배기//
  end when;

  //Thermodynamic Cycle Equations//
  if Period==61 then
    p=101.3e3;
    T=273+60;
  elseif Period==23 or Period==41 then
    der(m*u)=der(Q)-p*der(V)+h_in*der(m_in)+h_out*der(m_out);
    p*der(V)+V*der(p)=m*R_gas*der(T)+m*T*der(R_gas)+R_gas*T*der(m);
  elseif Period==12 then
    p=c*((thetta-thetta_1)^2)+d;
    T=e*((thetta-thetta_1)^2)+f;
  else
    p=c*((thetta-thetta_3)^2)+d;
    T=e*((thetta-thetta_3)^2)+f;
  end if;

  when Period==61 or Period==23 or Period==41 then
    i=0;
    p2=pre(p);
    T2=pre(T);
    c=0;
    d=0;
    e=0;
    f=0;
  elsewhen Period==12 then
    i=pre(i)+1;
    if i==1 then
      p2=(pre(p2))*(epsilon^k);
      T2=(pre(T2))*(epsilon^(k-1));
      c = (p2-pre(p2)) / (thetta_0-thetta_1);
      d = pre(p);
      e = (T2-pre(T2)) / (thetta_0-thetta_1);
      f = pre(T);
    else
      p2=pre(p2);
      T2=pre(T2);
      c=pre(c);
      d=pre(d);
      e=pre(e);
      f=pre(f);
    end if;
  elsewhen Period==34 then
    i=pre(i)+1;
    if i==1 then
      p2=(pre(p2))*(epsilon^(1-k));
      T2=(pre(T2))*(epsilon^(-k));
      c = (p2-pre(p2)) / (thetta_4-thetta_3);
      d = pre(p);
      e = (T2-pre(T2)) / (thetta_4-thetta_3);
      f = pre(T);
    else
      p2=pre(p2);
      T2=pre(T2);
      c=pre(c);
      d=pre(d);
      e=pre(e);
      f=pre(f);
    end if;
  end when;

  //Mass Conservation Law : m, m_in, m_out//
  m=m_in+m_out;
  m=m_fuel+m_air;
  m_air=m_fuel*AFR;

  if Period==61 then
    der(m_in)=m_dot_in;
    der(m_out)=0;
  elseif (Period==12 or Period==23 or Period==34) then
    der(m_in)=0;
    der(m_out)=0;
  else
    der(m_in)=0;
    der(m_out)=-m_dot_in;
  end if;

  //Cylinder Volume : V, V_c, V_s//
  V=V_c+V_s;
  V_c=V_s/(epsilon-1);
  V_s=(pi*(B^2)/4)*(L+R-R*cos(thetta_rad) + (( (L^2)-(R^2)*((sin(thetta_rad))^2)) ^(1/2)));
  S_p=(L+R-R*cos(thetta_rad)+(((L^2)-(R^2)*((sin(thetta_rad))^2))^(1/2)))
  annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
        coordinateSystem(preserveAspectRatio=false)));

  //Combustion sub-model : Q_x, Q_t//
  der(x_b)= if Period==23 then abs( 18000*Alpha*( (b+1)/del_thetta) *( ((thetta-thetta_0)/del_thetta)^b) *(1-x_b)) else 0;

  when thetta >= 720 then
    reinit(x_b, 0);
  end when;

  der(Q_x)=Q_t*der(x_b);
  Q_t=m_fuel*LHV_f;
  Alpha=log(1-x_star);

  //Thermodynamic model : Q_w, T, T_w//
  Q=Q_x-Q_w;

  Q_w_dot=der(Q_w);
  A_w=((pi*(B^2)/2)+(4*V/B))*(x_b^(1/2))+((pi*(B^2)/2)+(4*V/B))*0.001;
  A_w_max=(pi*(B^2)/2)+(4*V/B);
  Q_w_dot= Alpha_w*A_w_max*(T-T_w);

  WallHeatFlux = Alpha_w*T;

  //Woschni heat transfer coefficient//
  p_wos = p/(98*(10^3));
  p_c_wos = p_c/(98*(10^3));
  Alpha_w=110*( ((p_wos^0.8)*(w^0.8))/((T^0.53)*(B^0.2)));

  if (Period==23) or (Period==41) then
    w=2.28*C_m+3.24*(10^(-3))*((V_c*T)/(p*V))*(p_wos-p_c_wos); //Combustion, Expansion Period//
  elseif Period==61 then
    w=6.18*C_m; //Scavenging//
  else
    w=2.28*C_m; //Compression//
  end if;

  C_m=(S_p*omega)/30;
  p_c = 0.7*p;

  //Cooling Channel//
  Q_w_dot=U*A_w_max*T_LMTD;
  Q_w_dot=m_dot_c*Cpc*(T_co-T_ci);

  U=1/( (1/Alpha_w)+(1/h_c)+(t_w/k_w) );

  //Robust LMTD Method(T_LMTD varying the cases in the algorithm//
  del_T1 =  T - T_co;
  del_T2 =  T - T_ci;

  // 냉각수 측 상관식 //
  Nu_Dh=((h_c*D_h)/k_c);
  Nu_Dh=((f_c/8)*abs(Re_Dh-1000)*Pr)/(1+12.7*((f_c/8)^0.5)*((Pr^0.67)-1));
  D_h=(2*D_c*D_c)/(D_c+D_c);
  f_c=(0.790*log(Re_Dh)-1.64)^(-2);

  del_p_c=0.5*f_c*(L_cooling*((v_c^2)/(D_h*2*9.81)));

  Pr=(u_c*Cpc)/k_c;
  Re_Dh=((v_c) * D_h)/vis_c;
  vis_c=u_c/rho_c;
  v_c=m_dot_c/(rho_c*D_c*D_c);

 if (del_T1 > beta) and (del_T2 > beta) and abs(del_T1-del_T2)>tol then
    T_LMTD = (del_T1-del_T2) / ( (log(del_T1)) - (log(del_T2)));
  elseif (del_T1 > beta) and (del_T2 > beta) and abs(del_T1-del_T2)<tol then
    T_LMTD = (del_T1-del_T2) / 2;
  elseif (del_T1 > beta) and (del_T2 < beta) then
    T_LMTD = (del_T1 - beta) / ( (log(del_T1/beta)) * (1 - zeta * (del_T2 - beta)));
  elseif (del_T1 < beta) and (del_T2 > beta) then
    T_LMTD = (del_T2 - beta) / ( (log(del_T2/beta)) * (1 - zeta * (del_T1 - beta)));
  elseif (del_T1 < beta) and (del_T2 < beta) then
    T_LMTD = beta / ( (1 - zeta * (del_T1 - beta))  * (1 - zeta * (del_T2 - beta)));
  else
    T_LMTD = beta / ( (1 - zeta * (del_T1 - beta))  * (1 - zeta * (del_T2 - beta)));
 end if;

initial equation
  p = 101.3*(10^3);
  p2 = 101.3*(10^3);
  m_in= 100e-6;
  m=100e-6;
  T = 273+60;
  T2 = 273+60;
  thetta=0;
  x_b=0;
  c=0;
  d=0;
  e=0;
  f=0;
  i=0;

  annotation (experiment(StopTime=0.08));
end enginemodel_ContinuousEditEq;

1 Ответ

2 голосов
/ 20 марта 2020

Было бы неплохо иметь его в форме, которой могут управлять другие люди.

В любом случае, поскольку при t = 0 Period = 61, как насчет удаления всех других веток и просто протестировать модель? для 0 <тета <180? Это может помочь при отладке. </p>

И так как вы на нем: есть библиотека Testing, поставляемая с Dymola. Я рекомендую взглянуть на это. Идея может быть:

  • проверить модель без ветвления, для 0
  • сделать ссылку
  • реализовать ветвление для Period = 12 (180 < тета <340) </li>
  • тестовый период = 61 против эталона, а период = 12, исходя из ваших знаний
  • сделать ссылку на 0 <тета <340 </li>
  • реализовать ветвление для периода = ...
  • ... и так далее, пока вы не повторно внедрили (и не протестировали) все ветвление

Я знаю, что это не ответ на ваш вопрос, но я надеюсь это может помочь вам найти ответ самостоятельно. Это также работа, которая окупится в долгосрочной перспективе, когда вам придется вносить изменения в эту модель.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Похожие темы

...