1

Нужно написать программу, которая бы на выходе давала распределение температуры (T) по глубине(z), по выражению с интегралом (он должен быть вычислен численно), который я прикрепил картинкой. Не совсем понимаю каким образом это можно реализовать. Есть пример, но защитить его я вряд ли смогу, так что буду благодарен любой помощи: введите сюда описание изображения:

import math
import numpy as np
Q1 = 0.00058
Q2 = 0.00058
z1_к = 1000
z1_п = 1050
z2_к = 1500
z2_п = 1505
H_бнкт = 500
T_зж = 40
H_скв = 1505
H_уст = 0
T0 = 20
G = 0.03
c_ж = 4160
p_ж = 1000
Q = Q1+Q2
m_ж = 0.001004
r1 = 0.031
r2 = 0.0365
r3 = 0.057
r4 = 0.066
r5 = 0.108
c_гп = 880
tр = 86400
l_ж = 0.6
l_нкт = 43
l_ок = 43
l_ц = 0.73
l_мт = 0.6
l_гп = 0.8
dT0 = T_зж-T0
p_гп = 1500
Re = p_ж*Q/(np.pi*m_ж*(r1)**2)
Pr = m_ж*c_ж/l_ж
Nu = 0.021*Re**0.8*Pr**0.43
a_гп = l_гп/(c_гп*p_гп)
t = a_гп*tр/(r5)**2
Td = np.log(np.exp(-0.2*t)+(1.5-0.3719*np.exp(-t))*t**0.5)
h_НКТ = 1 
#h для ОК
h_ОК = 2 
#B для НКТ
B_НКТ = c_ж*p_ж*Q/h_НКТ
#B для ОК от БНКТ до 1 интервала перфорации
B_НКТ_1_инт = c_ж*p_ж*Q/h_ОК 
#B для ОК от 1 интервала перфорации до 2 интервала перфорации
B_1_инт_2_инт = c_ж*p_ж*Q2/h_ОК
N = 1000 #Количество разбиений по z от 0 до H_скв
h_z = (H_скв-H_уст)/N #Шаг по z
O_z = np.ones((0, N+1)) #Единичный массив  
Z = np.linspace(0,H_скв,num = N+1) 
T0 = 20
def B_глубина(z): #Определяю функцию B как, по сути, кусочную, тк прямой зависимости от z нет
    if 0 <= z <= H_бнкт:
        B = B_НКТ #постоянное числовое значение. Мне стоит сделать B массивом? 
    if H_бнкт < z <= z1_к:
        B = B_НКТ_1_инт #то же, что и в комментарии для B_нкт
    if z1_к < z <= z2_к:
        B = B_1_инт_2_инт #то же, что и в комментарии для B_нкт
    return B

def Integral():
    M=1000
    F=np.zeros(N+1)
    for i in range (1,N+1):
        Y=np.linspace (0,Z[i],num=M)
        for j in range(0,M):
            F[i]+=1/B_глубина(Y[j])
    return F
T = T0*np.exp(-Integral())
print (T)

Ошибка:

        runfile('C:/Users/Asus/Desktop/циклы.py', 
  wdir='C:/Users/Asus/Desktop')
  Traceback (most recent call last):

  File "C:\Users\Asus\Desktop\циклы.py", line 71, in <module>
  T = T0*np.exp(-Integral())

  File "C:\Users\Asus\Desktop\циклы.py", line 69, in Integral
    F[i]+=1/B_глубина(Y[j])

  File "C:\Users\Asus\Desktop\циклы.py", line 61, in B_глубина
    return B

UnboundLocalError: local variable 'B' referenced before assignment

Заранее спасибо за ответ!

0

1 ответ 1

1

Помогли разобраться - на всякий случай публикую получившийся код для тех, кому этот вопрос чем-то в теории может помочь:

import math
import numpy as np
import matplotlib.pyplot as plt
Q1 = 0.00058
Q2 = 0.00058
z1_к = 1000
z1_п = 1050
z2_к = 1500
z2_п = 1505
H_бнкт = 500
T_зж = 40
H_скв = 1505
H_уст = 0
T0 = 20
G = 0.03
c_ж = 4160
p_ж = 1000
Q = Q1+Q2
m_ж = 0.001004
r1 = 0.031
r2 = 0.0365
r3 = 0.057
r4 = 0.066
r5 = 0.108
c_гп = 880
tр = 86400
l_ж = 0.6
l_нкт = 43
l_ок = 43
l_ц = 0.73
l_мт = 0.6
l_гп = 0.8
dT0 = T_зж-T0
p_гп = 1500
Re = p_ж*Q/(np.pi*m_ж*(r1)**2)
Pr = m_ж*c_ж/l_ж
Nu = 0.021*Re**0.8*Pr**0.43
a_гп = l_гп/(c_гп*p_гп)
t = a_гп*tр/(r5)**2
Td = np.log(np.exp(-0.2*t)+(1.5-0.3719*np.exp(-t))*t**0.5)
#Р=h для НКТ
h_НКТ = 2*math.pi/(2/(l_ж*Nu) 
+1/l_нкт*np.log(r2/r1)+1/l_мт*np.log(r3/r2)+1/l_ок*np.log(r4/r3) 
+1/l_ц*np.log(r5/r4)+1/l_гп*Td)
#h для ОК
h_ОК = 2*math.pi/(2/(l_ж*Nu)+1/l_ок*np.log(r4/r3)+1/l_ц*np.log(r5/r4)+1/l_гп*Td)
#B для НКТ
B_НКТ = c_ж*p_ж*Q/h_НКТ
#B для ОК от БНКТ до 1 интервала перфорации
B_НКТ_1_инт = c_ж*p_ж*Q/h_ОК 
#B для ОК от 1 интервала перфорации до 2 интервала перфорации
B_1_инт_2_инт = c_ж*p_ж*Q2/h_ОК
N = 1000
Z = np.linspace(0, 1500, num = N+1) 
T0 = 20
M = 1000
def B_глубина(z): 
    if 0 <= z <= 500:
        B = B_НКТ 
    if 500 < z <= 1000:
        B = B_НКТ_1_инт
    if 1000 < z <= 1500:
        B = B_1_инт_2_инт
    return B

def Integral():
    F=np.zeros(N+1)
    for i in range (1,N+1):
        Y=np.linspace (0,Z[i],num=M)
        dz=Z[i]/M
        for j in range(0,M):
            F[i]+=dz/B_глубина(Y[j])
    return F
T = T0*np.exp(-Integral())
plt.plot(T,Z)
plt.gca().invert_yaxis()
plt.show()
print (T)

Ваш ответ

Нажимая «Отправить ответ», вы соглашаетесь с условиями пользования и подтверждаете, что прочитали политику конфиденциальности.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.