中国石油大学油藏数值模拟作业

刘月田老师的作业，大家可以参考参考。

import numpy as np
import math

xy = np.array（[[0 0 0 0 0 0 0 0 0 0 0]
               [0 259 222 200 190 180 185 0 0 0 0]
               [0 310 240 235 228 210 195 0 0 0 0]
               [0 330 290 270 250 230 205 0 0 0 0]
               [0 350 300 280 259 222 200 190 180 185 0]
               [0 340 320 290 310 240 235 228 210 195 0]
               [0 355 335 315 310 290 270 250 230 205 0]
               [0 0 0 0 325 300 280 240 210 215 0]
               [0 0 0 0 340 320 290 260 235 225 0]
               [0 0 0 0 355 335 315 295 275 255 0]
               [0 0 0 0 0 0 0 0 0 0 0]]）
# 将坐标换至以左下角为零点，横轴为x轴，纵轴为y轴的坐标系
xy = np.flip（xy 0）
xy = np.swapaxes（xy 0 1）


dt = 1  # h
dx = 200  # m
dy = 200  # m
C = 0.0002  # /MPa
VIS = 5  # mPa⋅s
RW = 10  # cm
RE = 0.208 * dx  # m
Q1 = 30  # m^3/d
pwf2 = 15  # MPa
pini = 20  # MPa
K = np.zeros（shape=（11 11））
H = np.zeros（shape=（11 11））
PRO = np.zeros（shape=（11 11））
for i in range（1 10）:
    for j in range（1 10）:
        if i < 7 and j > 3 or i > 3 and j < 7:
            K[i][j] = （xy[i][j] + 80） / 1000  # µm2
            H[i][j] = （xy[i][j] / 70） * 100  # cm
            PRO[i][j] = （xy[i][j] * 0.05 + 10） / 100
# 达西单位换算
dt = dt * 3600  # s
dx = dx * 100  # cm
dy = dy * 100  # cm
C = C * 10  # /10^-1MPa------------------------------------------------------------------------------------------------
VIS = VIS  # mPa⋅s
RW = RW  # cm
RE = RE * 100  # cm
Q1 = Q1 * 10 ** 6 / 86400  # cm^3/s
pwf2 = pwf2 * 10  # 10^-1MPa
pini = pini * 10  # 10^-1MPa

HKip = np.zeros（shape=（11 11））
HKis = np.zeros（shape=（11 11））
HKjp = np.zeros（shape=（11 11））
HKjs = np.zeros（shape=（11 11））
a = np.zeros（shape=（11 11））
b = np.zeros（shape=（11 11））
c = np.zeros（shape=（11 11））
d = np.zeros（shape=（11 11））
e = np.zeros（shape=（11 11））
qvr = np.zeros（shape=（11 11））
p = np.full（（11 11） float（pini））

for i in range（1 10）:
    for j in range（1 10）:
        if i < 7 and j > 3 or i > 3 and j < 7:
            HKip[i][j] = （2 * H[i + 1][j] * K[i + 1][j] * H[i][j] * K[i][j]） / （
                    H[i + 1][j] * K[i + 1][j] + H[i][j] * K[i][j]）
            HKis[i][j] = （2 * H[i - 1][j] * K[i - 1][j] * H[i][j] * K[i][j]） / （
                    H[i - 1][j] * K[i - 1][j] + H[i][j] * K[i][j]）
            HKjp[i][j] = （2 * H[i][j + 1] * K[i][j + 1] * H[i][j] * K[i][j]） / （
                    H[i][j + 1] * K[i][j + 1] + H[i][j] * K[i][j]）
            HKjs[i][j] = （2 * H[i][j - 1] * K[i][j - 1] * H[i][j] * K[i][j]） / （
                    H[i][j - 1] * K[i][j - 1] + H[i][j] * K[i][j]）

            a[i][j] =