如何从算法部分提取值？

Question

我正在使用抛物线PDE与代数方程混合，再加上所有这些方程。 我在OM（1.9.1）（“residualFcn [some number]”）中使用了Euler方法（Dassl太慢）和大容忍（用于快速模拟）和recive错误（对于这两种类型）。问题是求解器可以解决非线性系统（数学上，系统是正确的）。 第一个问题是什么类型的方法在OM中使用欧拉方法进行积分（显式或隐式或Crank-Nicholson..or ...）？所以我试图解决它的数值（显式欧拉方法（在代码下面的“新[N]”）（也许问题可以是CFL条件），但我有问题（样本重建的具体采样时间） 所以，第二个问题是指重现特定采样时间的值？！ 在下面的代码中有数组“a [3]”。 想法是针对每个“ts（采样时间）”重建3个节点中的值。 这怎么能 我怎样才能从算法部分当前值（即节点）传递到公式部分？ 我有.txt文件中的所有值，我用Matlab来绘制它们，但我不知道如何通过它，即在“方程”部分或其他方式 同样的问题是“new [N]”，针对特定的采样时间，节点的绘图函数（N）如果delta（t）/（delta（x））^ 2> = 0.5（delta（t）定义用户，并且参考公式部分，则delta（x）如下面的代码所示，并且空间离散化用于方程部分（经典前馈方法）），数值稳定性是否满足？ 同样的例子，但算法部分？ 问候

这里是代码：

model Euler1D 
import Modelica.Utilities.*; 
parameter Integer N=10; //50 
parameter Real Lp=1e-6; 
parameter Real deltax=1/(N-1)*Lp; 
Real a[3]; 
Real old[N]; 
Real new[N]; 
Real b; 
equation 
a[1]=if 
    (time>5) then 0 else time+5; 
a[2]=time; 
a[3]=2; 
when 
(sample(0,1)) then 
d=b; 
end when; 
algorithm 
// IN t=ts; 
when (sample(0,1)) then 
for i in 1:2 loop 
    b:=a[i]; 
Streams.print(String(time)+" "+String(a[i])+ "  "+String(b), "C:/Some_Path/text.txt"); 
end for; 
end when; 

// Another problem 
old[1]:=10; 
old [N]:=0; 
new[1]:=10; 
new [N]:=0; 
// Boundary 
for i in 2:N-1 loop 
old [i]:=10; 
new[i]:=10; 
end for; 

for dx in deltax:deltax:Lp-deltax loop // spatial discretization 

for i in 2:N-1 loop 
(new[i]):=(old[i]+0.5*(old[i + 1] +old[i-1]- 2*old[i])); 
//def:=def+abs(new[i]-old[i]); 
end for; 
for i in 2:N-1 loop 
old[i]:=new[i]; // switch the values 
end for; 

for i in 1:N loop 
Streams.print(String(time)+" "+ String(new[2]), "C:/Some_Path/Anel_Nodes.txt"); 
end for; 

annotation (uses(Modelica(version="3.2"))); 
end Euler1D; 

Answer 1

考虑您的意见，这是我会怎么构建模型：

model Euler1D 
    parameter Integer N=10 "Spatial discretization"; 
    parameter Modelica.SIunits.Length L=1; 
protected 
    parameter Modelica.SIunits.Length dx=L/(N-1) "Segment size"; 
    Real c[N] "Solution variable c"; 
    Real J[N] "Solution variable J"; 
    Real dc_dx "Spatial derivative at x==L"; 
    Real d2c_dx2[N] "Second spatial derivative of c"; 
initial equation 
    der(c[2:N-1]) = zeros(N-2); 
equation 
    // Equations for spatial derivatives 
    d2c_dx2[2:N-1] = { (c[i+1]-2*c[i]+c[i-1])/(dx^2) for i in 2:N-1}; 
    dc_dx = (c[N]-c[N-1])/dx; 

    // Boundary conditions 
    c[1] = 0 "Dirichlet B.C."; 
    dc_dx = 1 "Neumann B.C."; 

    // PDE 
    J = c .* c "J = c^2"; 
    der(c) = d2c_dx2+J "PDE"; 
end Euler1D; 

初始方程避免大的瞬态在模拟的开始。如此有效地，我提出的模型将为您提供稳定的状态解决方案。但是你可以使事物变化（例如边界条件）。

可能还存在一些错误，但我很难说，因为我对这个系统没有多少直觉，我也无法知道我使用的值是否是物理的有意义的。

但希望它概述了如何构建这样一个模型。这个模型运行在Dymola。我没有在OpenModelica中尝试过。所有这些都说了，你总是用PDE来解决这个问题，即时间积分方案不能明确看到物理学产生的隐式约束（例如CFL条件）。通常，可变时间步长算法具有误差估计，并且它们可以间接地看到由违反这些隐式约束而引入的不稳定性。但他们必须通过他们“摸索”。我不确定这将会成为上述模型的问题。

如果你想一点点的兴奋添加到您的解决方案，您可以更改公式

dc_dx = 1 "Neumann B.C."; 

...到...

dc_dx = 0.7+0.2*sin(time) "Neumann B.C."; 

...，你会得到一个不错的时间变化解决我注意到，在玩这个游戏的时候，如果你把梯度变得过于陡峭，模型会变得不稳定。无论这是基础数学方程的真正解决方案还是由于上述问题导致的数值假象，我都不知道。所以如果你想要一个稳定的解决方案，你可以放入这个模型是有限的。