ODE +用于热交换计算的偏方程

Question

我正在尝试编写代码t来解决流体1和2合并在一起时温度的微分方程，并考虑以下壁温：enter image description here

for dT1/dx  and dT2/dx i have decided to use Euler for mesch and they have been rewrite as:
(T1[i-1]-2*T1[i]+T1[i+1])/dx      and      (T2[i-1]-2*T2[i]+T2[i+1])/dx 

下面是定义数组的主要部分：

   T1=np.empty([n], dtype='float')
print('T1', T1)
print('T1shape', T1.shape)
print('T1dnim', T1.ndim)
T2=np.empty(([n]),dtype='float')
print('T2', T2)
print('T2shape', T2.shape)
print('T2dnim', T2.ndim)
Ts=np.empty([n],dtype='float')
print('Ts', Ts)
print('Tsshape', Ts.shape)
print('Tsdnim', Ts.ndim)enter code here

dTsdt = np.empty([n],dtype='float')
print('dTsdt', dTsdt)
print('dTsdtshape', dTsdt.shape)
print('dTsdtdnim', dTsdt.ndim)
dT1dt = np.empty([n],dtype='float')
print('dT1dt', dT1dt)
print('dT1dtshape', dT1dt.shape)
print('dT1dtdnim', dT1dt.ndim)
dT2dt = np.empty([n],dtype='float')
print('dT2dt', dT2dt)
print('dT2dtshape', dT2dt.shape)
print('dT2dtdnim', dT2dt.ndim)
t=np.linspace(1,10,n)
print('t', t)
print('tshape', t.shape)

并且集成已按如下方式完成：

for i in range(0, n):
    def mod(Ts, t):
        if i == 0:
            dTsdt[0] = ((T1[0] / tau1) - (T2[0] / tau2)) - (1 / tau2 - 1 / tau1) * Ts[0]
        elif i in range(1, n - 2):
            dTsdt[i] = ((T1[i,i] / tau1) - (T2[i,i] / tau2)) - (1 / tau2 - 1 / tau1) * Ts[i]
        else:

            dTsdt[n - 1] = ((T1[n-1,n - 1] / tau1) - (T2[n-1,n - 1] / tau2)) - (1 / tau2 - 1 / tau1) * Ts[n - 1]
        return dTsdt


    Ts0 = np.ones([n]) * 85
    Ts = odeint(mod, Ts0, dt)
    # Ts = Ts[:, -2]
    print('===============================================================')
    print('i', i)
    print('============================T surface =========================')
    print('Tss', Ts)
    print('Ts shape', Ts.shape)
#=================================================================================================
#================================ T1 Calculation =============================================
#------------------------------------------------------------------------------------------------
    def mod1(T1, t):
            if i == 0:
                dT1dt[0] = ((Ts[0,0]-T1[0])/tau1)-(u1/dx**2)*(T1s[0] - 2 * T1[0] + T1[0])
            elif i in range(1, n - 2):
                dT1dt[i] = ((Ts[i,i]-T1[i])/tau1)-(u1/dx**2) * (T1[i - 1]-2*T1[i]+T1[i + 1])
            else:
                dT1dt[n-1] = ((Ts[i,n-1]-T1[n-1])/tau1)-(u1/dx**2)*(T1[n-2]-2*T1[n-1]+T1s[n-1]
            return dT1dt


    t1 = np.linspace(0, 60, n)
    T10 = np.ones([n]) * 120
    T1 = odeint(mod1, T10, t1)
    # T1 = T1[:, -2]
    print('============================T1 Tempt===========================')
    print('T1', T1)
    print('T1type', T1.shape)
# ==================================================================================================
# ======================================= T2 Calculation ===========================================
# ---------------------------------------------------------------------------------------------------
    def mod2(T2, t):
            if i == 0:
                dT2dt[0] = ((Ts[0,0] - T2[0])/tau2 )- (u2/ dx**2) * (T1s[0]-2*T2[0]+T2[0])
            elif i in range(1, n - 2):
                dT2dt[i] = ((Ts[i,i] - T2[i])/tau2)-(u2/dx**2)*(T2[i-1]-2*T2[i] + T2[i+1])
            else:
                dT2dt[n-1] = ((Ts[i,n-1]-T2[n-1])/tau2)-(u2/dx**2)*(T2[n-2]-2*T2[n-1] +T1s[n-1])
            return dT1dt





    t1 = np.linspace(0, 1, n)
    T20 = np.ones([n]) * 140
    T2= odeint(mod2, T20, t1)
    #T2 = T2[:, -2]
    print('============================T2 Temp===========================')
    print('T2', T2)

您能不能帮助我纠正它，以某种方式不能给出正确的答案。

Answer 1

导数表达式似乎有错误：

(T1[i-1]-2*T1[i]+T1[i+1])/dx

这实际上是二阶导数乘以dx。一阶导数是

dT/dx = (T1[i+1] - T1[i-1])/(2*dx)

此外，对于偏微分方程，欧拉方法通常效果不佳；最佳方法取决于PDE的类型。

由于这是一个热流问题，因此您需要确保在每个离散的时间步长上都节省能量。我尚未验证您的特定实现是否如此，但通常，它有助于重写质量流量[kg / s]，比热[J /（kg-K）]和热量方面的方程式网格点之间的电导[W / K]。