円筒物体周りの流れによる強制空冷について ( その2 )¶
その1では、温度を仮定した1回の計算での算出について記載した. 温度が上昇する毎に、空気の物性値(動粘度、熱伝導度、密度、等)が変化するため、非線形的に解く必要がある.
物性値の変化¶
動粘度の温度依存性¶
-1.00000000e+02 5.83077000e-06
0.00000000e+00 1.35173000e-05
2.50000000e+01 1.58206000e-05
1.00000000e+02 2.35155000e-05
2.00000000e+02 3.54676000e-05
3.00000000e+02 4.91765000e-05
4.00000000e+02 6.45068000e-05
5.00000000e+02 8.13416000e-05
6.00000000e+02 9.96237000e-05
7.00000000e+02 1.19262000e-04
8.00000000e+02 1.40241000e-04
熱伝導度の温度依存性¶
-1.00000000e+02 1.57000000e-02
-5.00000000e+01 2.00000000e-02
-2.00000000e+01 2.24000000e-02
0.00000000e+00 2.41000000e-02
2.00000000e+01 2.57000000e-02
4.00000000e+01 2.72000000e-02
6.00000000e+01 2.87000000e-02
8.00000000e+01 3.02000000e-02
1.00000000e+02 3.16000000e-02
1.20000000e+02 3.31000000e-02
1.40000000e+02 3.45000000e-02
1.60000000e+02 3.59000000e-02
1.80000000e+02 3.72000000e-02
2.00000000e+02 3.86000000e-02
2.50000000e+02 4.18000000e-02
3.00000000e+02 4.49000000e-02
3.50000000e+02 4.79000000e-02
4.00000000e+02 5.08000000e-02
5.00000000e+02 5.62000000e-02
6.00000000e+02 6.13000000e-02
8.00000000e+02 7.09000000e-02
1.00000000e+03 8.02000000e-02
1.20000000e+03 8.91000000e-02
1.40000000e+03 9.70000000e-02
1.60000000e+03 1.04700000e-01
反復解法コード¶
コード¶
import os, sys
import numpy as np
import scipy as sp
# ========================================================= #
# === prepare__parameters === #
# ========================================================= #
def prepare__parameters( inpFile="dat/parameter.conf" ):
x_, y_ = 0, 1
# ------------------------------------------------- #
# --- [1] load parameter.conf --- #
# ------------------------------------------------- #
import nkUtilities.load__constants as lcn
params = lcn.load__constants( inpFile=inpFile )
# ------------------------------------------------- #
# --- [2] calculate dependent variables --- #
# ------------------------------------------------- #
params["pipe.area"] = params["pipe.diameter"]**2 * np.pi / 4
params["pipe.flow.m3_s"] = params["pipe.flow.L_min"] * 1e-3 / 60.0
params["fluid.velocity"] = params["pipe.flow.m3_s"] / params["pipe.area"]
params["target.area"] = 2.0*( np.pi / 4 * params["target.diameter"]**2 ) + np.pi * params["target.diameter"] * params["target.length"]
# ------------------------------------------------- #
# --- [3] load physical property --- #
# ------------------------------------------------- #
import nkUtilities.load__pointFile as lpf
tCondData = lpf.load__pointFile(inpFile=params["fluid.tConductivityFile"] )
viscoData = lpf.load__pointFile(inpFile=params["fluid.viscosityFile"] )
params["tConductivity_fit"] = sp.interpolate.interp1d( tCondData[:,x_], tCondData[:,y_], \
fill_value="extrapolate" )
params["viscosity_fit"] = sp.interpolate.interp1d( viscoData[:,x_], viscoData[:,y_], \
fill_value="extrapolate" )
# ------------------------------------------------- #
# --- [4] load HilpertModel Coefficient --- #
# ------------------------------------------------- #
coeffData = lpf.load__pointFile(inpFile=params["fluid.HilpertCoeffFile"] )
min_, max_, c_, m_ = 0, 1, 2, 3
eps = 1.e-8
coeffData[:,max_] = coeffData[:,max_] - eps
xD = np.concatenate( [ coeffData[:,min_][None,:], coeffData[:,max_][None,:] ], axis=0 )
cD = np.concatenate( [ coeffData[:,c_ ][None,:], coeffData[:,c_ ][None,:] ], axis=0 )
mD = np.concatenate( [ coeffData[:,m_ ][None,:], coeffData[:,m_ ][None,:] ], axis=0 )
xD = np.reshape( xD, (-1,) )
cD = np.reshape( cD, (-1,) )
mD = np.reshape( mD, (-1,) )
params["HilpertCoeff_c_fit"] = sp.interpolate.interp1d( xD, cD, fill_value="extrapolate" )
params["HilpertCoeff_m_fit"] = sp.interpolate.interp1d( xD, mD, fill_value="extrapolate" )
return( params )
# ========================================================= #
# === display__variables === #
# ========================================================= #
def display__variables( params=None, skip_keys=None ):
# ------------------------------------------------- #
# --- [1] arguments check --- #
# ------------------------------------------------- #
if ( params is None ): sys.exit( "[display__variables.py] params == ???" )
if ( skip_keys is None ):
skip_keys = [ "control.iterMax", "control.maxResidual", "control.verbose", \
"tConductivity_fit", "viscosity_fit", \
"HilpertCoeff_c_fit", "HilpertCoeff_m_fit" ]
keys = list( set( params.keys() ) - set( skip_keys ) )
# ------------------------------------------------- #
# --- [2] display --- #
# ------------------------------------------------- #
print( "\n" + "-"*70 )
print( " {0:>30} :: {1:>30} ".format( "key", "value" ) )
print( "-"*70 )
for key in keys:
print( " {0:>30} :: {1:>30} ".format( key, params[key] ) )
print( "-"*70 + "\n" )
return()
# ========================================================= #
# === main.py === #
# ========================================================= #
def main( parameterFile=None ):
# ------------------------------------------------- #
# --- [1] arguments check & load parameters --- #
# ------------------------------------------------- #
if ( parameterFile is None ): parameterFile="dat/parameter.conf"
params = prepare__parameters( inpFile=parameterFile )
print( "\n" + "="*90 )
print( "[solver__forcedCoolingAroundCylinder.py] Begin Main Loop" )
print( "="*90 + "\n" )
# ------------------------------------------------- #
# --- [2] Main Loop --- #
# ------------------------------------------------- #
for ik in range( params["control.iterMax"] ):
# ------------------------------------------------- #
# --- [4-1] update physical property --- #
# ------------------------------------------------- #
params["fluid.thermalConductivity"] = params["tConductivity_fit"]( params["target.T"] )
params["fluid.viscosity"] = params["viscosity_fit"] ( params["target.T"] )
# ------------------------------------------------- #
# --- [4-2] update cooled temperature --- #
# ------------------------------------------------- #
params["target.Told"] = params["target.T"]
params["fluid.Re"] = params["fluid.velocity"] * params["target.diameter"] / params["fluid.viscosity"]
params["coeff.c"] = params["HilpertCoeff_c_fit"]( params["fluid.Re"] )
params["coeff.m"] = params["HilpertCoeff_m_fit"]( params["fluid.Re"] )
params["fluid.Nu"] = params["coeff.c"] * params["fluid.Re"]**params["coeff.m"] * params["fluid.prandtl"] ** ( 1./3. )
params["heat_transfer"] = params["fluid.Nu"] * params["fluid.thermalConductivity"] / params["target.diameter"]
params["target.dT"] = params["target.heatload"] / ( params["heat_transfer"] * params["target.area"] )
params["target.T"] = params["fluid.temperature_infty"] + params["target.dT"]
params["target.T"] = params["target.Told"] + params["control.relaxation"] * ( params["target.T"] - params["target.Told"] )
# ------------------------------------------------- #
# --- [4-3] check convergence --- #
# ------------------------------------------------- #
residual = np.abs( params["target.T"] - params["target.Told"] ) / ( params["target.Told"] )
if ( residual < params["control.maxResidual"] ):
converged = True
print( "\n" + "="*90 )
print( " Reach convergence. at #. of iteration :: {}".format( ik ) )
print( "="*90 + "\n" )
print( " ( iteration, temperature, residual ) == ( {0:8}, {1:10.4f}, {2:10.4e} )".format( ik, params["target.T"], residual ) )
display__variables( params=params )
print( "\n" + "="*90 + "\n" )
break
# ------------------------------------------------- #
# --- [4-4] display variables --- #
# ------------------------------------------------- #
if ( params["control.verbose"] ):
print( "\n" + "-"*80 )
print( " iteration :: {}".format( ik ) )
print( "-"*80 + "\n" )
print( " ( iteration, temperature, residual ) == ( {0:8}, {1:10.4f}, {2:10.4f} )".format( ik, params["target.T"], residual ) )
display__variables( params=params )
print( "\n" + "-"*80 + "\n" )
else:
print( " ( iteration, temperature, residual ) == ( {0:8}, {1:10.4f}, {2:10.4f} )".format( ik, params["target.T"], residual ) )
if ( not( converged ) ):
print( "\n" + "[solver__forcedCoolingAroundCylinder.py] does not converged... [CAUTION] " + "\n" )
# ========================================================= #
# === Execution of Pragram === #
# ========================================================= #
if ( __name__=="__main__" ):
parameterFile = "dat/parameter.conf"
main( parameterFile=parameterFile )
パラメータファイル¶
# for single code ( formula )
fluid.viscosity float 4.9e-5 # m2/s
fluid.thermalConductivity float 4.5e-2 # @ 300 degree
fluid.prandtl float 0.7
fluid.temperature_infty float 20.0 # K or degree,
pipe.diameter float 12e-3 # m
pipe.flow.L_min float 70 # L/min
target.diameter float 13e-3 # m
target.length float 100e-3 # m
target.heatload float 100 # W
target.T float 110 # K or degree
coeff.c float 0.193
coeff.m float 0.618
control.iterMax integer 100
control.maxResidual float 1.e-5
control.verbose logical False
control.relaxation float 0.5
fluid.tConductivityFile string dat/thermalConductivity.dat
fluid.viscosityFile string dat/viscosity.dat
fluid.HilpertCoeffFile string dat/HilpertModel.dat
実行結果¶
Reference¶
[1] Website of BUILDING PHYSICS RESEARCH GROUP, "空気の物性値(温度による)" ( https://lee-lab.net/blog-contents-037/ )