LBM学习记录5 Python实现IB二维圆柱绕流 1.0

浸润边界（IB）下二维圆柱绕流，bug1.0版

import numpy as np
from math import sqrt
from numpy import pi as PI
from numpy import sin,cos
import matplotlib.pyplot as plt
from matplotlib import cm
import matplotlib
matplotlib.use('Agg')

# 读取拉格朗日点
def read_IB_Points(struct_name):filename = struct_name+'.ibpoint'with open(filename) as f:Nb = int(f.readline().strip())xLag,yLag = np.loadtxt(f,unpack=True) # unpack将每一列当成一个向量输出return(Nb,xLag,yLag)

######### def ##############

# 找近邻邻居点
def give_1D_NoneDelta_Indices(xy_Lag,N_xy,d_xy):delta = 4ind = np.floor(xy_Lag/d_xy + 1) # np.floor向下取整indices = np.tile(ind,(delta,1)).T # np.tile扩展数组indices += -delta/2+1+np.arange(delta)indices = (indices-1)%N_xyreturn indices

def give_NonZero_Delta_Indices_XY(xLag,yLag,Nx,Ny,dx,dy):delta = 4xInds = give_1D_NoneDelta_Indices(xLag,Nx,dx)yInds = give_1D_NoneDelta_Indices(yLag,Ny,dy)xInds = np.tile(xInds,(1,delta))yInds = np.tile(yInds,(1,delta))return(xInds.astype('int'),yInds.astype('int'))

def give_Eulerian_Lagrangian_distance(x,y,L):distance = np.abs(x-y)distance = np.minimum(distance,L-distance)return distance

def give_Delta_Kernel(x,dx):RMAT = np.abs(x)/dxdelta = RMATrow,col = x.shapefor ii in range(row):for jj in range(col):r = RMAT[ii,jj]if r<1:delta[ii,jj] = ( (3 - 2*r + sqrt(1 + 4*r - 4*r*r) ) / (8*dx) )elif (r<2) and (r>=1):delta[ii,jj] = ( (5 - 2*r - sqrt(-7 + 12*r - 4*r*r) ) / (8*dx) )else:delta[ii,jj] = 0return delta

def give_Me_Perturbed_Distance(U,V,dx,dy,delta_X,delta_Y,xInds,yInds):mat_X = U[xInds,yInds]*delta_X*delta_Ymat_Y = V[xInds,yInds]*delta_X*delta_Ymove_X = mat_X.sum(1) * (dx*dy)move_Y = mat_Y.sum(1) * (dx*dy)return(move_X, move_Y)

# step 1 更新拉格朗日边界位置
def please_Move_Lagrangian_Point_Positions(x,y,Nx,Ny,Lx,Ly,dx,dy,Nb,ds,xLag,yLag,U,V,xL_P,yL_P,dt):delta = 4Nx = NxNy = NyLx = LxLy = Lydx = dxdy = dyNb = Nbds = dsU = UV = Vdt = dtxInds, yInds = give_NonZero_Delta_Indices_XY(xLag,yLag,Nx,Ny,dx,dy)xLH_aux = xLag % LxyLH_aux = yLag % LyxL_H_ReSize = np.tile(xLH_aux,(delta**2,1)).TyL_H_ReSize = np.tile(yLH_aux,(delta**2,1)).T#     if(np.isscalar(xL_P)):distX = give_Eulerian_Lagrangian_distance(x[xInds],xL_H_ReSize,Lx)distY = give_Eulerian_Lagrangian_distance(y[yInds],yL_H_ReSize,Ly)
#     else:
#         distX = give_Eulerian_Lagrangian_distance(x[xInds],xL_H_ReSize,Lx)
#         distY = give_Eulerian_Lagrangian_distance(y[yInds],yL_H_ReSize,Ly)delta_X = give_Delta_Kernel(distX,dx)delta_Y = give_Delta_Kernel(distY,dy)move_X, move_Y = give_Me_Perturbed_Distance(U,V,dx,dy,delta_X,delta_Y,xInds,yInds)xL_Next = xL_P + dt * move_XyL_Next = yL_P + dt * move_Yreturn(xL_Next,yL_Next)

########### def #####

def give_Me_Delta_function_approximations_For_Force_Calc(Nx,Ny,Lx,Ly,dx,dy,Nb,x,y,xLag,yLag):delta = 4Nx = NxNy = NyLx = LxLy = Lydx = dxdy = dyNb = NbindX = give_1D_NoneDelta_Indices(xLag, Nx, dx)indY = give_1D_NoneDelta_Indices(xLag, Ny, dy)indLagAux = np.arange(Nb)ind_Lag = np.tile(indLagAux,(delta,1)).T.astype('int')distX = give_Eulerian_Lagrangian_distance(x[indX.astype('int')],xLag[ind_Lag],Lx)distY = give_Eulerian_Lagrangian_distance(y[indY.astype('int')],yLag[ind_Lag],Ly)# 初始化delta_X和Y矩阵delta_X = np.zeros((Nb,Nx))delta_Y = np.zeros((Ny,Nb))delta_1D_x = give_Delta_Kernel(distX,dx)delta_1D_y = give_Delta_Kernel(distY,dy)delta_X[ind_Lag,indX.astype('int')] = delta_1D_xdelta_Y[indY.astype('int'),ind_Lag] = delta_1D_yreturn(delta_X,delta_Y)

# step 2 计算来自拉格朗日点的力
def please_Find_Lagrangian_Forces_On_Eulerian_grid(Nx,Ny,Lx,Ly,dx,dy,Nb,ds,xLag,yLag,x,y):delta = 4Nx = NxNy = NyLx = LxLy = Lydx = dxdy = dyNb = Nbds = dsfx = np.zeros(Nb)fy = np.zeros(Nb)# 保存拉格朗日力F_Lag = np.empty((fx.size,2))F_Lag[:,0] = fxF_Lag[:,1] = fydelta_X,delta_Y = give_Me_Delta_function_approximations_For_Force_Calc(Nx,Ny,Lx,Ly,dx,dy,Nb,x,y,xLag,yLag)fxds = np.diag(fx*ds)fyds = np.diag(fy*ds)Fx = delta_Y @ fxds @ delta_XFy = delta_Y @ fyds @ delta_Xreturn(Fx.T,Fy.T,F_Lag)

############ def ########

LBM

def give_Me_Macroscopic_Rho_U(fin, Fx, Fy, dt,Nx,Ny,v):rho = np.sum(fin, axis=0)u = np.zeros((2,Nx,Ny))for i in range(9):u[0,:,:] = (v[i,0] * fin[i,:,:])/rho + Fx*dt/(2*rho)u[1,:,:] = (v[i,1] * fin[i,:,:])/rho + Fy*dt/(2*rho)U = u[0,:,:]V = u[1,:,:]return(rho,u,U,V)

# 平衡态计算
def equilibrium(rho, u, Nx, Ny,v,t):usqr = 3/2 * (u[0]**2 + u[1]**2)feq = np.zeros((9, Nx, Ny))for i in range(9):cu = 3 * (v[i,0]*u[0,:,:] + v[i,1]*u[1,:,:])feq[i,:,:] = rho*t[i] * (1 + cu + 0.5*cu**2 - usqr)return feq

# 初始速度曲线：几乎为零，有一个轻微的扰动来触发不稳定。
def inivel(d, x, y):return (1-d) * 0.04 * (1+1e-4*sin(y/2*2*PI))

def please_Exe_LBM(Nx,Ny,fin,vel,Fx,Fy,dt,v,t,omega):col1 = np.array([0, 1, 2])col2 = np.array([3, 4, 5])col3 = np.array([6, 7, 8])fin[col3, -1, :] = fin[col3, -2, :]rho,u,U,V = give_Me_Macroscopic_Rho_U(fin,Fx,Fy,dt,Nx,Ny,v)u[:, 0, :] = vel[:, 0, :]U[0,:] = vel[0,0,:]V[0,:] = vel[1,0,:]rho[0,:] = 1/(1-u[0,0,:]) * (np.sum(fin[col2,0,:], axis=0) + 2*np.sum(fin[col3,0,:], axis=0))feq = equilibrium(rho,u,Nx,Ny,v,t)fin[[0,1,2],0,:] = feq[[0,1,2],0,:]+fin[[8,7,6],0,:]-feq[[8,7,6],0,:]fout = fin - omega * (fin - feq)for i in range(9):fin[i, :, :] = np.roll(np.roll(fout[i,:,:],v[i,0], axis=0), v[i,1], axis=1)return(fin,vel,U,V)

def main():struct_name = 'cylinder'Time = 3.5dt = 0.005current_time = 0i = 1Nx = 40Ny = 20Lx = 40Ly = 20dx = Lx/Nxdy = Ly/Nyds = Lx/(2.*Nx)cx, cy, r = Nx//4, Ny//2, Ny//9  # 圆柱坐标Re = 220.0         # 雷诺数uLB = 0.04        # 模型的入口速度nulb = uLB*r/Re   # 粘性系数omega = 1/(3*nulb+0.5)  # 松弛频率ly = Ny-1         # 用于计算入口，为模型增加扰动，当雷诺数较小时，计算缺少扰动v = np.array([[1,1], [1,0], [1,-1],[0,1], [0,0], [0,-1],[-1,1], [-1,0], [-1,-1]])t = np.array([1/36, 1/9, 1/36,1/9, 4/9, 1/9,1/36, 1/9, 1/36])# 创建欧拉格点x = np.arange(0, Lx, dx)y = np.arange(0, Ly, dy)X,Y = np.meshgrid(x, y)# 读取拉格朗日节点Nb,xLag,yLag = read_IB_Points(struct_name)xLag_P = xLag #初始化前拉格朗日点yLag_P = yLagU = np.zeros((Nx,Ny)) # 初始化速度V = np.zeros((Nx,Ny))vel = np.fromfunction(inivel, (2, Nx, Ny))fin = equilibrium(1,vel,Nx,Ny,v,t)while current_time < Time:# step1:更新拉格朗日点位置xLag_h,yLag_h = please_Move_Lagrangian_Point_Positions(x,y,Nx,Ny,Lx,Ly,dx,dy,Nb,ds,xLag,yLag,U,V,xLag,yLag,dt)
#         print('xLag_h')
#         print(xLag_h)# step2:计算来自拉格朗日点的力Fx,Fy,F_Lag = please_Find_Lagrangian_Forces_On_Eulerian_grid(Nx,Ny,Lx,Ly,dx,dy,Nb,ds,xLag,yLag,x,y)
#         print('Fx')
#         print(Fx)# step3:LBMfin,vel,U,V = please_Exe_LBM(Nx,Ny,fin, vel, Fx, Fy, dt, v, t, omega)
#         print('U')
#         print(U[15])#step4:更新拉格朗日点的位置xLag_P = np.array(xLag_h)yLag_P = np.array(yLag_h)xLag,yLag = please_Move_Lagrangian_Point_Positions(x,y,Nx,Ny,Lx,Ly,dx,dy,Nb,ds,xLag,yLag,U,V,xLag_P,yLag_P,dt)
#         print('xLag_P')
#         print(xLag_P)
#         print('xLag')
#         print(xLag)current_time += dt#         wait = input('Press enter to contiune...')#可视化xlabel = str(current_time)+' s'plt.figure(figsize=(8, 4))
#         plt.clf()plt.xlabel(xlabel)plt.axis('equal')plt.plot(xLag,yLag,'r-',xLag,yLag,'*')plt.imshow(np.sqrt(U**2+V**2).transpose(), cmap=cm.Reds)plt.savefig("/share/home/yszhang/PBS/LBM/IB_D2Q9/pic/"+str(i)+".jpg")i += 1

main()

import cv2file_dir = '/share/home/yszhang/PBS/LBM/IB_D2Q9/pic/'
video = cv2.VideoWriter('video.avi',cv2.VideoWriter_fourcc(*'MJPG'),5,(1280,720))for i in range(1,701):img = cv2.imread(file_dir+str(i)+'.jpg')     img = cv2.resize(img,(1280,720))video.write(img)video.release()

LBM学习记录5 Python实现IB二维圆柱绕流 1.0相关推荐

fluent二维叶型仿真_ICEM划分嵌套网格之二维圆柱绕流
首先,介绍一下嵌套网格.网络上关于嵌套网格的的内容大多数是关于直接利用软件进行计算的过程,而对于前处理过程中的网格生成过程并没有什么描述,其实这种技术已经在学术界流传已久,只是用的都是自己的程序算法, ...
matlab圆柱饶流,有限元法解二维圆柱绕流问题.pdf
目录目录 1 1．问题描述 1 2．相关的有限元理论基础 1 2.1 二次泛函极值原理和里兹解法 1 2.2 伽辽金加权余数法 2 3．二维圆柱绕流的有限元解法 ...
LBM学习记录3 Python实现D2Q9圆柱绕流
在https://blog.csdn.net/weixin_58913471/article/details/117561995?spm=1001.2101.3001.6650.1&utm_m ...
java学习记录十五：集合二Collections、Set、Map
java学习记录十五:集合二一.Collections工具类一.解释二.常用方法 1.打乱集合顺序 2.按照默认规则排序 3.按指定规则排序 4.批量添加元素二.可变参数一.解释二.写法 ...
Python生成动态二维码，运用神库：qrcode
一.介绍 1.1 二维码二维码又称二维条码,常见的二维码为 QR Code,QR 全称 Quick Response.是一个近几年来移动设备上超流行的一种编码方式,在现在的生活中二维码随处可见.我们 ...
python myqr制作二维码生成器_用Python生成动态二维码，只要5行代码，拥有你的个性二维码！...
原标题:用Python生成动态二维码,只要5行代码,拥有你的个性二维码! 前言文的文字及图片来源于网络,仅供学习.交流使用,不具有任何商业用途,版权归原作者所有,如有问题请及时联系我们以作处理. P ...
如何用python制作动态二维码，提升表白成功率？
来源:凹凸数据本文约1000字,建议阅读5分钟. 本文教你用python制作动态二维码,助你表白成功! 关注数据派THU(DatapiTHU)后台回复"20200520"获取完整 ...
深入浅出python机器学习_如何用python画(绘制)二维函数(二维图)?
参考文档 python 如何绘制二维函数? from matplotlib import pyplot as plt import numpy as np low=lambda x:10000 if ...
Python+OpenCV：二维直方图(2D Histograms)
Python+OpenCV:二维直方图(2D Histograms) ################################################################# ...

LBM学习记录5 Python实现IB二维圆柱绕流 1.0

LBM学习记录5 Python实现IB二维圆柱绕流 1.0相关推荐

最新文章

热门文章