管道压降及fluent 第一层网格高度yplus y+计算

# -*- coding: utf-8 -*-
"""
Created on Fri Jan 20 11:24:22 2023

@author: PS
摩擦阻力计算式：ΔP1 = λ×(L/D)×(ρv^2)/2
用来计算管道的沿程压降,湍流，雷诺数大于4000
一般定义圆管内雷诺数大于4000为湍流，2000-4000为过度流，小于2000为层流
"""

from scipy.optimize import fsolve
from math import sqrt,log10
import numpy as np
import math

#输入参数
ccd = 2e-3  #绝对粗糙度m
L =15 #管道长度m
D = 200e-3  #管道直径m,特征尺寸，4*A/P,A为面积m3,P为周长
den = 998.2 #介质密度kg/m3 水998.2kg/m3  空气：1.225kg/m3
ve = 3 #管道截面的速度为m/s
dn = 0.001003   #dn 动力黏性系数η(Pa·s)  空气：1.7894e-05 水：0.001003
yplus = 0.8 #期望的y+
r = 1.2 #网格增长率
# n = 30 #网格层数


#公式输入e:当量粗糙度e(绝对粗糙度(m))、管道直径D(m)、雷诺数，求解阻力摩擦系数
# Colebrook-White方程
def f1_solve(e,D,Re):
    sol = fsolve(lambda x:
                 1/sqrt(x)+2*log10(e/3.7/D+2.51/(Re*sqrt(x)))
                 , np.array([0.001]))
    print('Colebrook-White方程阻力摩擦系数f为：',f"{round(sol[0],6)}")
    return sol[0]

# Swamee-Jain 方程
def f2_solve(e,D,Re):
    f = 0.25/pow(log10(e/D/3.7+5.74/pow(Re,0.9)),2)
    print('Swamee-Jain方程阻力摩擦系数f为：',f"{round(f,6)}")
    return f

# 实用流体阻力手册
def f3_solve(Re):
    f = 1/pow(1.8*log10(Re)-1.64,2)
    print('实用流体阻力手册沿程阻力摩擦系数（光滑圆管）f为：',f"{round(f,6)}")
    return f



#  密度den(kg/m^3),速度vel(m/s),直径d(m),dn 动力黏性系数η(Pa·s)，求解雷诺数
def re(den,vel,d,dn):
    res = den*vel*d/dn
    print('当前计算的雷诺数Re为：',f"{round(res,3)}")
    return res
  
#求压降    
Re = re(den,ve,D,dn) #雷诺数
pa = f1_solve(ccd, D, Re)*(L/D)*den*pow(ve,2)/2  #压降-表压
pa2 = f2_solve(ccd, D, Re)*(L/D)*den*pow(ve,2)/2  #压降-表压
print("-------------------------")
print('Colebrook-White方程计算的压降为：',f"{round(pa,5)} Pa")
print('Swamee-Jain方程计算的压降为：',f"{round(pa2,5)} Pa")
print("-------------------------")

#湍流强度,通常I<1%称为低湍流强度,大于10%称为高湍流强度，5%左右通常称为中等湍流强度。
I = 0.16*pow(Re,-1/8)
print('湍流强度：',f"{round(I*100,3)} %")

#湍动能 公式：https://baike.baidu.com/item/%E6%B9%8D%E6%B5%81%E5%8B%95%E8%83%BD/4317324
k = 3/2*pow(ve*I,2)
print('湍动能 ：',f"{round(k*100,3)}% m2/s2")


#求解壁面切应力
Cf= pow(1/(-3.6*np.log10(6.9/Re+pow(ccd/3.7/D,1.11))),2)#壁面摩擦系数 skin coff 
tw = 0.5*den*pow(ve,2)*Cf #壁面切应力 wall shear stress
ut = pow(tw/den,0.5) #壁面摩擦速度 m/s

#求解yplus
ydis = yplus*dn/(ut*den) #第一层网格距离m,应该乘以2，因为体积有限法是计算的网格中心

if Re <5e5:
    h_layer = 4.91*D/pow(Re,0.5)
else:
    h_layer = 0.38*D/pow(Re,1/5)
    
# 以5为底数 print(math.log(25, 5))	# 2.0
ds = 1-h_layer*(1-r)/ydis
n = math.log(ds,r) #最少网格层数
yt = ydis*(1-pow(r,n))/(1-r) #边界层网格的总厚度
yfin = ydis * pow(r,n-1)  #最外层网格厚度
print("-------------------------")
print('壁面摩擦系数：',f"{round(Cf,6)}")
print('壁面切应力：',f"{round(tw,5)} Pa")
print("-------------------------")
print('第一层网格距离：',f"{round(ydis*2,10)} m")
print('最少网格层数：',f"{round(n,0)}")
print('边界层厚度：',f"{round(h_layer,10)} m")
print('边界层网格的总厚度：',f"{round(yt,10)} m")
print('最外层网格厚度：',f"{round(yfin,10)} m")

当前计算的雷诺数Re为： 597128.614
Colebrook-White方程阻力摩擦系数f为： 0.038006
Swamee-Jain方程阻力摩擦系数f为： 0.038076
-------------------------
Colebrook-White方程计算的压降为： 17071.81727 Pa
Swamee-Jain方程计算的压降为： 17103.15089 Pa
-------------------------
湍流强度： 3.035 %
湍动能 ： 1.243% m2/s2
-------------------------
壁面摩擦系数： 0.015668
壁面切应力： 70.37927 Pa
-------------------------
第一层网格距离： 7.5683e-06 m
最少网格层数： 31.0
边界层厚度： 0.0053161841 m
边界层网格的总厚度： 0.0053161841 m
最外层网格厚度： 0.0008891842 m

计算摩擦系数f参考如下图片：