Метод конечных разностей

1. Подготовка

Вводим физические параметры материала



In [1]:

    
# Коэффициент теплопроводности
lam = 401
# Коэффициент удельной теплоемкости
c = 385
# Плотность материала
ro = 8900

# Вычисляем коэффициент задающий физические свойства материала
alpha = lam/(c*ro)

Теперь параметры задачи



In [2]:

    
# Длина стержня в метрах
stick_length = 1
# Время исследования в секундах
time_span = 200

# Шаг по длине стержня
dx = 0.02
# Шаг по времени
dt = 5

Вводим граничные и начальные условия



In [3]:

    
# Начальная температура стержня
Tinit = 300
# Температура справа
Tright = 350
# Температура слева
Tleft = 320

Введем дополнительные переменные для удобства



In [4]:

    
# Вычислим количество шагов по времени
time_step_count = int(time_span / dt) + 1
# Вычислим количество шагов по координате
len_step_count = int(stick_length / dx) + 1

# Размерность гломальной матрицы
dim_global_matrix = time_step_count*len_step_count

Определим дополнительную функцию, для преобразования локальных координат, в индекс в глобальной матрицы



In [5]:

    
def get_gidx(time_line, x):
    return time_line*len_step_count+x

2. Решение

на этой стадии мы сформируем матрицы, применим шаблон, и решим систему уравнений



In [6]:

    
# Подключим математическую библиотеку
import numpy as np

Введем коэффициенты для шаблона

Значения уравнений были получены для уравнения Пуассона, с применением двухслойного шеститочечного шаблона, в зависимости от значения sigma будут 3 различных шаблона:

0 - Явная схема (Не забывайте о условиях устойчивости)
0.5 - Схема Кранка-Николсона
1 - Не явная четырехточечная схема



In [7]:

    
sigma = 1
c = 2*alpha*dt*(1-sigma) - dx**2 # u p m
v = 2*dt * alpha*sigma + dx**2 # u p+1 m
h = -alpha*dt*sigma # u p+1 m-1
j = -alpha*dt*sigma # u p+1 m+1
k = -alpha*dt*(1- sigma) # u p m-1
p = -alpha*dt*(1- sigma) # u p m+1



In [8]:

    
# Формируем матрицы
A = np.diag(np.ones(dim_global_matrix))
B = np.zeros([dim_global_matrix])

Заполним вектор B, введем начальные и граничные условия



In [9]:

    
# Начальное условие
B[0:len_step_count]=Tinit

# Граничные условия
B[len_step_count::len_step_count] = Tleft
right_column = len_step_count+len_step_count-1
B[right_column::len_step_count] = Tright

Формируем глобальную матрицу, заполняем нашим шаблоном



In [10]:

    
for time_line in range(1, time_step_count):
    for x_index in range(1, len_step_count-1):
        active_row = get_gidx(time_line, x_index)
        A[active_row, get_gidx(time_line, x_index)] = v
        A[active_row, get_gidx(time_line, x_index+1)] = h
        A[active_row, get_gidx(time_line, x_index-1)] = j
        A[active_row, get_gidx(time_line-1, x_index)] = c
        A[active_row, get_gidx(time_line-1, x_index-1)] = k
        A[active_row, get_gidx(time_line-1, x_index+1)] = p

Решаем СЛАУ



In [11]:

    
X =  np.linalg.solve(A, B)



In [12]:

    
print(X)









    



[ 300.          300.          300.         ...,  342.5959792   346.28151105
  350.        ]

3. Результаты

Построим 2D графики для каждого временного слоя



In [13]:

    
%matplotlib inline

import matplotlib.pyplot as plt

x_axis_step = np.arange(0, stick_length+dx/2.0, dx)
for n in range(1, time_step_count):
    plt.plot(x_axis_step, X[len_step_count*n:len_step_count*n+len_step_count])

plt.show()

Строим 3D график



In [22]:

    
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

fig = plt.figure()
ax = fig.gca(projection='3d')

XA = np.arange(0, stick_length+dx/2.0, dx)
YA = np.arange(0, time_span+dt/2.0, dt)
XA, YA = np.meshgrid(XA, YA)
ZA = X.reshape((time_step_count,len_step_count)).tolist()

surf = ax.plot_surface(XA, YA, ZA, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()
#mpld3.show()



In [23]:

    
# Показать интерактивно 
%matplotlib qt 
from matplotlib import interactive
interactive(True)

fig = plt.figure()
ax = fig.gca(projection='3d')

XA = np.arange(0, stick_length+dx/2.0, dx)
YA = np.arange(0, time_span+dt/2.0, dt)
XA, YA = np.meshgrid(XA, YA)
ZA = X.reshape((time_step_count,len_step_count)).tolist()

surf = ax.plot_surface(XA, YA, ZA, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)
fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()



In [25]:

    
### Сохранить результат в файл



In [26]:

    
f=open('results.csv','wb')
header = ";".join([ 'X='+ str(x) for x in x_axis_step ])
f.write(  bytes(header+'\n', 'UTF-8'))
np.savetxt(f,ZA, delimiter=";")
f.close()