1

我基本上是在尝试使用 GEKKO 求解动态二阶偏微分方程。这是供参考的方程： 2-D 传热方程。

这里，t 是时间，T 是温度，(k, rho, Cp,le, sigma 和 Z) 都是常数。

T 已作为数组输入（我认为这就是问题所在）。

为了对方程进行数值求解，需要输入其他方程的边界条件。

这是代码：

import numpy as np
from gekko import GEKKO

m = GEKKO()

tf = 5*60*60
dt = int(tf/300)+1
m.time = np.linspace(0,tf,dt)

# Number of nodes
n = 41

# Length of domain
Lx = 1
Ly = Lx # square domain

x_div = np.linspace(0,Lx,n)
y_div = np.linspace(Ly,0,n)

[X, Y] = np.meshgrid(x_div, y_div)

# step size
dx = x_div[1] - x_div[0]
dy = y_div[1] - y_div[0]

# Temp. initialization
T = m.Var(np.ones((n,n))*290)

# Equation set-up

# Middle segments
for i in range(1,n-1):
    for j in range(1,n-1):
        m.Equations(T[i,j].dt() == (k/(rho*Cp)*((T[i+1,j]-2*T[i,j]+T[i-1,j])/dx**2 + (T[i,j+1]-2*T[i,j]+T[i,j-1])/dy**2))
                    + (G_total - em*sigma*(T[i,j]**4-T_surr**4))/(rho*Cp*Z))

# Boundary Conditions
m.Equations(T[0,:]==310,
           T[-1,:]==310,
           T[1:-2,0]==315,
           T[1:-2,-1]==315,
           
           T[0,0]==312,
           T[n,0]==312,
           T[0,n]==312,
           T[n,n]==312)

基本上，我正在尝试解决这个由温度组成的网格。我收到以下错误：“numpy.float64”对象没有属性“dt”

如果我只写 T 而不是 T[i,j]，我会得到这个错误：'int' object is not subscriptable

我的问题：

  1. GEKKO 是否能够求解这些本质上是二维的方程？我该怎么做？
  2. 为此目的，还有其他很酷的库吗？随着时间的推移，我需要能够绘制具有板温度的等高线图；为此我需要解方程。

感谢您的时间。

4

1 回答 1

1

这是 Lutz Lehmann 建议用于m.Array定义的解决方案T。

热解

import numpy as np
from gekko import GEKKO

m = GEKKO(remote=False)

tf = 5*60*60
dt = int(tf/300)+1
#m.time = np.linspace(0,tf,dt)
# for testing
m.time = [0,0.01,0.02,0.03]

# Number of nodes
n = 41

# Length of domain
Lx = 1
Ly = Lx # square domain

# Define for material
k = 1; rho = 8000; Cp = 500
G_total=1; em=1; sigma=1
T_surr=298; Z=1

x_div = np.linspace(0,Lx,n)
y_div = np.linspace(Ly,0,n)

[X, Y] = np.meshgrid(x_div, y_div)

# step size
dx = x_div[1] - x_div[0]
dy = y_div[1] - y_div[0]

# Temp. initialization
T = m.Array(m.Var,(n,n),value=290)

# Equation set-up

# Middle segments
for i in range(1,n-1):
    for j in range(1,n-1):
        m.Equation(rho*Cp*T[i,j].dt() == (k*\
                                    ((T[i+1,j]-2*T[i,j]+T[i-1,j])/dx**2 \
                                     + (T[i,j+1]-2*T[i,j]+T[i,j-1])/dy**2))
                    + (G_total - em*sigma*(T[i,j]**4-T_surr**4))/Z)

# Boundary Conditions
m.Equations([T[0,i]==310 for i in range(1,n-1)])
m.Equations([T[-1,i]==310 for i in range(1,n-1)])
m.Equations([T[i,0]==315 for i in range(1,n-1)])
m.Equations([T[i,-1]==315 for i in range(1,n-1)])
m.Equations([T[0,0]==312, T[n-1,0]==312, \
             T[0,n-1]==312, T[n-1,n-1]==312])

m.options.IMODE=7
m.solve(disp=False)

import matplotlib.pyplot as plt

plt.figure(figsize=(8,8))
for i in range(0,4):
    for j in range(0,4):
        plt.subplot(4,4,i*4+j+1)
        plt.plot(m.time,T[i,j].value)
plt.savefig('heat.png',dpi=600)
plt.show()

还有其他使用 Gekko 的双曲线和抛物线 PDE示例。Gekko 的优势在于优化。对于模拟，最好使用 PDE 的标准模拟方法。此外，您可以使用IMODE=7.

于 2021-07-05T01:38:01.663 回答