mathmodel/dev1.py at main · shoyim/mathmodel · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma

def filtration():
    k0 = 1e-13
    myu0 = 1e-2
    P0 = 5e+5
    bet_y = 3e-10

    tMax = 3600

    lv = 1000
    lp = 500
    betta = 1
    alf = 0.7

    bet = 1 + betta

    tau = 100
    h = 0.05
    L = 40
    T = int(tMax / tau) + 1
    N = int(L / h) + 1

    # Initial and boundary conditions
    P = np.zeros((T, N))
    P[:, 0] = P0
    P[0, :] = 0
    P[:, -1] = 0
    U = np.zeros((T, N))
    kap = k0 / (myu0 * bet_y)

    for i in range(N):
        P[1, i] = P[0, i]

    kv = lv / ((tau ** bet) * gamma(3 - bet))
    kp = lp / ((tau ** alf) * (h ** 2) * gamma(2 - alf))

    kv1 = lv / ((tau ** betta) * gamma(2 - betta))  # velocity
    kp2 = lp / ((tau ** alf) * gamma(2 - alf))  # velocity

    A = kap / h ** 2 + kap * kp
    B = 1 / tau + kv + 2 * kap / h ** 2 + 2 * kap * kp
    C = kap / h ** 2 + kap * kp
    alpha = np.zeros(N)
    beta = np.zeros(N)
    beta[1] = P0

    for j in range(1, T - 1):
        for i in range(1, N - 1):
            sv = 0
            for k in range(1, j):
                sv += (P[k + 1, i] - 2 * P[k, i] + P[k - 1, i]) * ((j - k + 1) ** (2 - bet) - (j - k) ** (2 - bet))

            sp1 = 0
            sp2 = 0
            sp3 = 0
            st = 0
            for k in range(j):
                sp1 += (P[k + 1, i + 1] - P[k, i + 1]) * ((j - k + 1) ** (1 - alf) - (j - k) ** (1 - alf))
                sp2 += (P[k + 1, i] - P[k, i]) * ((j - k + 1) ** (1 - alf) - (j - k) ** (1 - alf))
                sp3 += (P[k + 1, i - 1] - P[k, i - 1]) * ((j - k + 1) ** (1 - alf) - (j - k) ** (1 - alf))
                st += (U[k + 1, i] - U[k, i]) * ((j - k + 1) ** (1 - betta) - (j - k) ** (1 - betta))  # velocity

        for i in range(1, N - 1):
            F = (P[j, i] / tau - kv * sv + 2 * kv * P[j, i] - kv * P[j - 1, i] + kp * kap * sp1 - kap * kp * P[j, i + 1] - 2 * kap * kp * sp2 + 2 * kap * kp * P[j, i] + kap * kp * sp3 - kap * kp * P[j, i - 1])
            alpha[i + 1] = C / (B - A * alpha[i])
            beta[i + 1] = (F + A * beta[i]) / (B - A * alpha[i])

        P[j + 1, -1] = 0
        for i in range(N - 2, -1, -1):
            P[j + 1, i] = alpha[i + 1] * P[j + 1, i + 1] + beta[i + 1]

        for i in range(N - 1):
            U[j + 1, i] = abs((-k0 * (P[j + 1, i + 1] - P[j + 1, i] + kp2 * (sp1 + P[j + 1, i + 1] - P[j, i + 1] - sp2 - P[j + 1, i] + P[j, i])) / (myu0 * h * (1 + kv1)) - (kv1 * st - kv1 * U[j, i]) / (1 + kv1)))

        U[j + 1, -1] = U[j + 1, -2]  # velocity

    # Plotting
    plt.plot(np.arange(N) * h, P[T - 1, :])
    plt.show()

filtration()