运用龙格库塔法解大雷洛数平板绕流问题

# -*- coding: utf-8 -*-
"""
Created on Mon Dec 16 19:40:52 2019

@author: Franz
"""

# 将高阶微分方程降为常微分方程组
# F'''+0.5F''F=0
# 方程组
# F= x
# F'= X' = y
# F'' = Y' = z
# F''' = Z' = -0.5F''F=-0.5yx


# define the calculate the space = 10000

import numpy as np
import matplotlib.pyplot as plt


plt.rcParams['figure.figsize'] = (6.0, 4.0)
plt.rcParams['figure.dpi'] = 200 #分辨率

h = 0.01;
error_less = 1e-6;
max_step = 10000;

x = np.zeros(max_step+1);
y = np.zeros(max_step+1);
z = np.zeros(max_step+1);

x[0] = 0;
y[0] = 0;
z[0] = 1;

K = np.zeros(4);
L= np.zeros(4);
M = np.zeros(4);

'''
f = y
g = z
z = -0.5xz
'''
for i in range(max_step):
    K[0] = h*y[i];   
    L[0] = h*z[i];
    M[0] = -0.5*h*x[i]*z[i];
    
    K[1] = h*(y[i]+0.5*L[0]);    
    L[1] = h*(z[i]+0.5*M[0]);
    M[1] = h*(-0.5*(x[i]+0.5*K[0])*(z[i]+0.5*M[0]));
    
    K[2] = h*(y[i]+0.5*L[1]);
    L[2] = h*(z[i]+0.5*M[1]);
    M[2] = h*(-0.5*(x[i]+0.5*K[1])*(z[i]+0.5*M[1]));
    
    K[3] = h*(y[i]+0.5*L[2]);
    L[3] = h*(z[i]+0.5*M[2]);
    M[3] = h*(-0.5*(x[i]+K[2])*(z[i]+M[2]));
    
    x[i+1] = x[i] + (K[0]+2*K[1]+2*K[2]+K[3])/6;
    y[i+1] = y[i] + (L[0]+2*L[1]+2*L[2]+L[3])/6;
    z[i+1] = z[i] + (M[0]+2*M[1]+2*M[2]+M[3])/6;
    
    A = y[i+1]**(-1.5);
    if (y[i+1]-y[i]) <= error_less:
        break;

x1 = np.zeros(max_step+1);
y1 = np.zeros(max_step+1);
z1 = np.zeros(max_step+1);
x1[0] = 0;
y1[0] = 0;
z1[0] = A;

index = 0;
for i in range(max_step):
    K[0] = h*y1[i];   
    L[0] = h*z1[i];
    M[0] = -0.5*h*x1[i]*z1[i];
    
    K[1] = h*(y1[i]+0.5*L[0]);    
    L[1] = h*(z1[i]+0.5*M[0]);
    M[1] = h*(-0.5*(x1[i]+0.5*K[0])*(z1[i]+0.5*M[0]));
    
    K[2] = h*(y1[i]+0.5*L[1]);
    L[2] = h*(z1[i]+0.5*M[1]);
    M[2] = h*(-0.5*(x1[i]+0.5*K[1])*(z1[i]+0.5*M[1]));
    
    K[3] = h*(y1[i]+0.5*L[2]);
    L[3] = h*(z1[i]+0.5*M[2]);
    M[3] = h*(-0.5*(x1[i]+K[2])*(z1[i]+M[2]));
    
    x1[i+1] = x1[i] + (K[0]+2*K[1]+2*K[2]+K[3])/6;
    y1[i+1] = y1[i] + (L[0]+2*L[1]+2*L[2]+L[3])/6;
    z1[i+1] = z1[i] + (M[0]+2*M[1]+2*M[2]+M[3])/6;
    
    index = i+1;
    if (y1[i+1]-y1[i]) <= error_less:
        break;

'''
定义物性:
    Ue = 5 m/s
    ratio = miu / rho = 1.49e-5
'''

Ue = 5;
ratio = 1.49e-5;

# 转化宏观参量
y0 = np.zeros((50,index+1));
u = np.zeros((50,index+1));
v = np.zeros((50,index+1));
x0 = np.zeros((50,index+1));
for i in np.arange(0,50):
    for j in range(0,index+1):
        x0[i,j] = i*0.01;
        g1=np.sqrt(x0[i,j]*2*ratio/Ue);
        g2=np.sqrt(ratio*Ue/2/(x0[i,j]+1e-6));
        eta = j*h;
        y0[i,j] = g1*eta;
        u[i,j]=Ue*y1[j];
        v[i,j]=g2*(eta*y1[j]-x1[j]); 


plt.figure();
for i in [20,30,40]:
    plt.plot(y0[i,:],u[i,:],label = 'x = %2.2f'%(x0[i,j]));

plt.legend(loc='upper left');
plt.ylabel('u');
plt.xlabel('y');
plt.savefig('./y_vs_u.png')
plt.show();

plt.figure();
for i in [20,30,40]:
    plt.plot(y0[i,:],v[i,:],label = 'x = %2.2f'%(x0[i,j]));
plt.legend(loc='upper left');
plt.ylabel('v');
plt.xlabel('y');
plt.savefig('./y_vs_v.png')
plt.show();

M = np.hypot(u, v);
plt.figure();
a = plt.contourf(x0,y0,M,50,origin='lower',cmap='jet');
plt.colorbar(a, ticks=np.max(M)*np.arange(0,1.1,0.1))
plt.savefig('./contour_velocity.png')
plt.show();


with open('./plate_flow.txt','w') as f:
    f.write('x	y	u	v
');
    for i in range(50):
        for j in range(index+1):
            f.write(str(x0[i,j])+'	'+str(y0[i,j])+'	'+str(u[i,j])+'	'+
                    str(v[i,j])+'
');