python - 使用 GEKKO 求解二维热方程
问题描述
我基本上是在尝试使用 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
我的问题:
- GEKKO 是否能够求解这些本质上是二维的方程?我该怎么做?
- 为此目的,还有其他很酷的库吗?随着时间的推移,我需要能够绘制具有板温度的等高线图;为此我需要解方程。
感谢您的时间。
解决方案
这是 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
.
推荐阅读
- ssl - 是否有解组 ASN.1 格式的 TLS 证书的接口
- ios - 在 Swift 中使用 API Gateway AWSTask 对象的正确方法是什么
- python - 在其他列的条目中制作列表长度的熊猫数据框列
- python - 在同一级别但也从顶层导入脚本
- components - 迭代时更改lwc组件中的表达式
- javascript - 为什么我得到控制台错误 findIndex() is not a function
- javascript - 隐藏按钮 12423
- version - Yocto CI 内部版本号?PR 服务不增加 ${PR}
- android - 错误:运行子进程cordova时发生错误
- php - 不使用 sql 外键的缺点