Реализация 1D тепловой диффузии PDE в Modelica (Dymola)

Я пытался реализовать пример одномерной диффузии тепла из «Объектно-ориентированного моделирования и моделирования с помощью Modelica 3.3» Питера Фрицона на основе подхода конечной разности функций сетки к одномерной диффузии тепла, но я получаю сообщение об ошибке:

Translation of Heat_diffusion_Test_1D.HeatDiffusion1D:

The DAE has 50 scalar unknowns and 50 scalar equations.

Cannot find differentiation function:
Heat_diffusion_Test_1D.DifferentialOperators1D.pder(
u, 
1)
with respect to time

    Failed to differentiate the equation
    Heat_diffusion_Test_1D.FieldDomainOperators1D.right(Heat_diffusion_Test_1D.DifferentialOperators1D.pder (
    u, 
    1)) = 0;

    in order to reduce the DAE index.

    Failed to reduce the DAE index.

    Translation aborted.

    ERROR: 1 error was found

    WARNING: 1 warning was issued

Я совсем новый пользователь языка программирования, у вас есть идеи, как я могу решить эту проблему? Видимо проблема в краевом условии правой границы. Заранее спасибо!

вот код:

package Heat_diffusion_Test_1D
import SI = Modelica.SIunits;
  model HeatDiffusion1D
    import SI = Modelica.SIunits;
    parameter SI.Length L=1;
    import Heat_diffusion_Test_1D.Fields.*;
    import Heat_diffusion_Test_1D.Domains.*;
    import Heat_diffusion_Test_1D.FieldDomainOperators1D.*;
    import Heat_diffusion_Test_1D.DifferentialOperators1D.*;
    constant Real PI=Modelica.Constants.pi;

    Domain1D rod(left=0, right=L, n=50);
    Field1D u(domain=rod, redeclare type FieldValueType=SI.Temperature, start={20*sin(PI/2*x)+ 300 for x in rod.x});

  equation 
    interior(der(u.val))=interior(4*pder(u,x=2));
    left(u.val)=20*sin(PI/12*time)+300;
   right(pder(u,x=1)) = 0;
    //right(u.val)=320;

    annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
          coordinateSystem(preserveAspectRatio=false)));
  end HeatDiffusion1D;

  package Fields

    record Field1D "Field over a 1D spatial domain"
      replaceable type FieldValueType = Real;
      parameter Domains.Domain1D domain;
      parameter FieldValueType start[domain.n]=zeros(domain.n);
      FieldValueType val[domain.n](start=start);

      annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
            coordinateSystem(preserveAspectRatio=false)));
    end Field1D;
    annotation ();
  end Fields;

  package Domains
    import SI = Modelica.SIunits;
    record Domain1D "1D spatial domain"
      parameter SI.Length left=0.0;
      parameter SI.Length right=1.0;
      parameter Integer n=100;
      parameter SI.Length dx = (right-left)/(n-1);
      parameter SI.Length x[n]=linspace(right,left,n);

      annotation (Icon(coordinateSystem(preserveAspectRatio=false)), Diagram(
            coordinateSystem(preserveAspectRatio=false)));
    end Domain1D;
    annotation ();
  end Domains;

  package FieldDomainOperators1D 
    "Selection operators of field on 1D domain"

    function left "Returns the left value of the field vector v1"
      input Real[:] v1;
      output Real v2;

    algorithm 
      v2 := v1[1];
    end left;

    function right "Returns the left value of the field vector v1"
      input Real[:] v1;
      output Real v2;

    algorithm 
      v2 := v1[end];

    end right;

    function interior 
      "returns the interior values of the field as a vector"
      input Real v1[:];
      output Real v2[size(v1,1)-2]; 

    algorithm 
      v2:= v1[2:end-1];

    end interior;


  end FieldDomainOperators1D;

  package DifferentialOperators1D 
    "Finite difference differential operators"

    function pder 
      "returns vector of spatial derivative values of a 1D field"
      input Heat_diffusion_Test_1D.Fields.Field1D f;
      input Integer x=1 "Diff order - first or second order derivative";
      output Real df[size(f.val,1)];

    algorithm 
      df:= if x == 1 then SecondOrder.diff1(f.val, f.domain.dx)
     else 
         SecondOrder.diff2(f.val, f.domain.dx);

    end pder;

    package SecondOrder 
      "Second order polynomial derivative approximations"

      function diff1 "First derivative"
        input Real u[:];
        input Real dx;
        output Real du[size(u,1)];

      algorithm 
        du := cat( 1, {(-3*u[1] + 4*u[2] - u[3])/2/dx},  (u[3:end] - u[1:end-2])/2/dx, {(3*u[end] -4*u[end-1] + u[end-2])/2/dx});

      end diff1;

      function diff2
        input Real u[:];
        input Real dx;
        output Real du2[size(u,1)];
      algorithm 
        du2 :=cat( 1, {(2*u[1] - 5*u[2] + 4*u[3]- u[4])/dx/dx}, (u[3:end] - 2*u[2:end-1] + u[1:end-2])/dx/dx, {(2*u[end] -5*u[end-1] + 4*u[end-2] - u[end-3])/dx/dx});
      end diff2;
      annotation ();
    end SecondOrder;

    package FirstOrder 
      "First order polynomial derivative approximations"

      function diff1 "First derivative"
        input Real u[:];
        input Real dx;
        output Real du[size(u,1)];

      algorithm 
                // Left, central and right differences
        du := cat(1, {(u[2] - u[1])/dx}, (u[3:end]-u[1:end-2])/2/dx, {(u[end] - u[end-1])/dx});

      end diff1;
      annotation ();
    end FirstOrder;
    annotation ();
  end DifferentialOperators1D;
  annotation (uses(Modelica(version="3.2.1")));
end Heat_diffusion_Test_1D;

person John_3265    schedule 11.05.2016    source источник
comment
Я не уверен, что оператор pder определен в текущей версии стандартной библиотеки modelica, и я не думаю, что библиотека 3.3, на которую ссылается книга, была выпущена?   -  person Christoph    schedule 11.05.2016
comment
@Christoph, pder не определен в текущей версии, это правильно. Я определяю это как функцию в коде. Стандартная библиотека Modelica — 3.2.1. Спецификация языка Modelica версии 3.3 уже реализована в Dymola (2016 FD01).   -  person John_3265    schedule 11.05.2016


Ответы (1)


Ваша реализация преобразуется и моделируется в Dymola 2015 FD01, если вы inline выполняете все функции в пакете FieldDomainOperators1D и DifferentialOperators1D с annotation(Inline=true). Причина указана Питером Фрицсоном на странице 406:

[...] Обратите внимание, что семантика inline необходима для функции pder() и функции сетки. [...] Это позволяет размещать вызовы такой функции в местах, где вызовы функций обычно не допускаются. [...]

person Lukas Exel    schedule 12.05.2016
comment
Я не уверен, правильно ли я понимаю этот ответ. Мне было интересно, если бы вы были настолько ребенком, чтобы добавить правильную версию кода в свой пост? - person Foad; 06.08.2017