LB 学习日记

LB学习日记

---

圆柱绕流算法

# -*- coding: utf-8 -*-
"""
Created on Mon Dec  9 12:46:49 2019

@author: Franz
"""

import numpy as np
import matplotlib.pyplot as plt
import numba as nb

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

# function
# 1st:equilibrum function
@nb.jit(nopython=True)
def Feq(k,rho,u):
    eu = e[k,0]*u[0] + e[k,1]*u[1];
    uv = u[0]**2 + u[1]**2;
    feq = rho*w[k]*(1.0+3.0*eu+4.5*eu*eu-1.5*uv);
    return feq

# 2nd:evolution: including move and collide,calculate micro-parameter
@nb.njit()
def Evolution(f,rho,u):
    F = np.zeros((NX+1,NY+1,Q));
    for i in range(0,NX+1):
        for j in range(1,NY):
            for k in range(Q):
                ip = i - e[k,0];
                jp = j - e[k,1];
                if (ip>=0) and (ip<=NX):
                    F[i,j,k] = f[ip,jp,k] + (Feq(k,rho[ip,jp],u[ip,jp])-
                     f[ip,jp,k])/tau_f;

    u0 = np.empty((NX+1,NY+1,2));
    for i in range(1,NX):
        for j in range(1,NY):
            u0[i,j,0] = u[i,j,0];
            u0[i,j,1] = u[i,j,1];
            rho[i,j] = 0;
            u[i,j,0] = 0;
            u[i,j,1] = 0;
            for k in range(Q):
                f[i,j,k] = F[i,j,k];
                rho[i,j] += f[i,j,k];
                u[i,j,0] += e[k,0]*f[i,j,k];
                u[i,j,1] += e[k,1]*f[i,j,k];
            u[i,j,0] /= (rho[i,j]+1e-30);
            u[i,j,1] /= (rho[i,j]+1e-30);
    # left & right marging
    for j in range(1,NY):
        u[0,j,0] = U*(1.0+1e-4*np.sin(j/NY*2*np.pi));
        rho[0,j] = rho[1,j];
        for k in range(Q):
            f[0,j,k]=Feq(k,rho[0,j],u[0,j])+f[1,j,k]-Feq(k,rho[1,j],u[1,j]);

    for j in range(1,NY):
        for k in range(Q):
            rho[NX,j] = rho[NX-1,j];
            u[NX,j] = u[NX-1,j];
            f[NX,j,k] = Feq(k,rho[NX,j],u[NX,j])+f[NX-1,j,k]-Feq(k,rho[NX-1,j],u[NX-1,j]);
    # top & bottom margin
    for i in range(NX+1):
        for k in range(Q):
            if (k == 4) or (k ==7) or (k==8):
                f[i,0,k] = f[i,NY,k]
            if (k == 2) or (k ==5) or (k==6):
                f[i,NY,k] = f[i,0,k]

    # for cylinder margin
    for i in range(0,NX):
        for j in range(0,NY):
            if index[i,j] == True:
                u[i,j,0] = 0;
                u[i,j,1] = 0;
            sumnum = 0;
            for k in range(Q):
                ip = i - e[k,0];
                jp = j - e[k,1];
                if index[ip,jp] == True:
                    sumnum += 1;
            if sumnum < 9 & sumnum > 1 & index[i,j]==True:
                if i >= cx & j >= cy:
                    rho[i,j] = rho[i+1,j+1];
                    f[i,j,k] = Feq(k,rho[i,j],u[i,j])+f[i+1,j+1,k]-Feq(k,rho[i+1,j+1],
                     u[i+1,j+1]);
                if i >= cx & j<=cy:
                    rho[i,j] = rho[i+1,j-1];
                    f[i,j,k] = Feq(k,rho[i,j],u[i,j])+f[i+1,j-1,k]-Feq(k,rho[i+1,j-1],
                     u[i+1,j-1]);
                if i < cx & j < cy:
                    rho[i,j] = rho[i-1,j-1];
                    f[i,j,k] = Feq(k,rho[i,j],u[i,j])+f[i-1,j-1,k]-Feq(k,rho[i-1,j-1],
                     u[i-1,j-1]);
                if i < cx & j >=cy:
                    rho[i,j] = rho[i-1,j+1];
                    f[i,j,k] = Feq(k,rho[i,j],u[i,j])+f[i-1,j+1,k]-Feq(k,rho[i-1,j+1],
                     u[i-1,j+1])
                
    return f,u,u0


# model parameter
global Q,NX,NY,U;
Q = 9;
NX = 480;
NY = 160;
U = 0.1;

global e,w,tau_f,index;
e = np.array([[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]],dtype=int);
w = np.array([4.0/9,1.0/9,1.0/9,1.0/9,1.0/9,1.0/36,1.0/36,1.0/36,1.0/36]);


# main
# init
dx = 1.0;
dy = 1.0;
Lx = dx * NX;
Ly = dy * NY;
dt = dx;
c = dx/dt;
rho0 = 1.0;
Re = 150;
cx = 50; # the (x,y)-cylinder center
cy = 85;
R = 15;
niu = U * 2 * R / Re;
tau_f = 3.0*niu + 0.5;
# cylinder
index = np.zeros((NX+1,NY+1),dtype=bool)
for i in range(NX+1):
    for j in range(NY+1):
        if (i-cx)**2+(j-cy)**2 <= R**2:
            index[i,j] = True;

rho = np.zeros((NX+1,NY+1),dtype='float')
u = np.zeros((NX+1,NY+1,2),dtype='float');

f = np.zeros((NX+1,NY+1,Q));

for i in range(NX+1):
    for j in range(NY+1):
        u[0,j,0] = U*(1.0+1.0e-4*np.sin(j/NY*2*np.pi)) ;
        rho[i,j]=rho0;
        for k in range(Q):
            f[i,j,k] = Feq(k,rho[i,j],u[i,j]);# 圆柱部分也进行了平衡分布函数的计算，且其部分不为零

max_steps = 8000;
for i in range(max_steps+1):
    f,u,u0 = Evolution(f,rho,u);
    ux = u[:,:,0];
    uy = u[:,:,1];
    x,y = np.meshgrid(np.arange(np.size(ux,0)),np.arange(np.size(ux,1)));
    print('The run times is %d, u|nx/2 = %f'%(i,ux[NX//2,NY//2]));
    M = np.hypot(ux, uy)
    if i % 50 == 0:
        plt.figure();
        a = plt.quiver(ux[::12,::12].T,uy[::12,::12].T,M[::12,::12])
        plt.figure();
        levels = np.linspace(0,np.max(M),100)
        plt.contourf(M.T,levels,origin='lower',cmap='jet')
        plt.savefig('./_cylinder_'+str(i)+'.png')
        plt.show()