はじめに
研究やアルバイトを通して、レーザー光を物体にあててその熱分布をシュミレーションしたいといったことがあったので計算を行いました。
計算量の問題で一次元の熱伝導の計算となっています。
伝熱方程式
一次元の熱伝導方程式は下記のように表せます。
$$ \frac{\partial T}{\partial t} = a \frac{\partial^2 T}{\partial x^2}$$
ここで、\(T,t,x,a\)はそれぞれ、温度\((K)\)、時間\((s)\)、長さ\((m)\)、熱拡散率\((m^2/s)\)を表しています。
これを計算機を用いて、計算するためには有限要素に分けて計算する必要があります。このため中心差分法を用いて離散化すると以下のように書き直すことができます。
$$ \frac{T^{n+1}(x_i)-T^n (x_i)}{\Delta t}= a \frac{T^n (x_i + \Delta x) + T^n (x_i – \Delta x) – 2T^n (x_i)}{\Delta x^2} $$
これを整理すると
$$ T^{n+1}(x_i) = d (T^n (x_i + \Delta x) + T^n (x_i – \Delta x)) + (1-2d) T^n (x_i) $$
\(d = a \Delta t /\Delta x^2\)であり、安定条件のために以下の関係を満たす必要があります。
$$ d < \frac{1}{2}$$
$$ \Delta t < \frac{\Delta x^2}{2a}$$
つまり、空間の分解能と時間の分解能はトレードオフの関係にあり、計算機の性能に制限がある限りどちらの分解能を求めるか考える必要があります。
熱の拡散の様子
まず、初めに一定の熱量がデルタ関数的に入射してきたときのことを考えます。
今回の計算では\( Sc_2 O_3\)といった結晶にレーザー光を入射させています。
今回の仮定では
- 結晶中で画像のようにガウス分布を持ったレーザー光が入射される
- 結晶端面は十分に冷却されている
- 結晶長は2 mm、レーザーの直径は50 µm
として計算を行いました。コードは以下のようになります。
"""
module import
"""
import matplotlib.pyplot as plt
import numpy as np
import plotly.graph_objects as go
"""
変数定義
"""
r0 = 25 #um -beam radius
xmin,xmax = 0,2000 #um - crystal length
t = 1 #sec - max time
Nt = 10**5 #sample point(time)
Nx = 10**2 #sample point(length)
kp = 16.5/0.683/4*10**6 #熱拡散係数(um^2/s)← 熱拡散係数(m^2/s) = 熱伝導率(W/mK)/(比熱 (J/(g K)) 密度(g/m^3))
TC = 16.5*10**-6 #熱伝導率(W/mK)
dx = (xmax-xmin)/Nx #1step time
dt = (t)/Nt #1step length
Q = 1.5 #power
abs = 0.5 #absorption rate
rep_rate = 1000 #Hz -reptation rate
x = np.linspace(xmin,xmax,Nx) #空間の行列
Tnew = np.zeros(Nx) #初期の温度
print("満たさなければいけない条件(2kpdt<dx**2)",2*kp*dt,"<",dx**2)
"""
数式定義
"""
#ガウス分布
def Gaussian(x,center,r0):
return np.exp(-2*((x-center)/(r0))**2)
#熱伝導方程式
def T(Tx_0,Tx_1,Tx_2,kp,dx,dt):
return Tx_1 + kp*dt*(Tx_0 -2*Tx_1 +Tx_2)/(dx**2)
"""
熱拡散の計算
"""
#最初に熱が入射
Tbefor = Gaussian(x,x[int(Nx/2)],r0)
#温度の分布情報と時間を保存しておく
cen = []
ti = []
#ある時間になったら画像を書いてもらいたいためその時間の指定
fre = [0,2,20,200,2000,20000,200000]
for j in range(Nt):
#print(Tnew)
for i in range(Nx-2):
x_0, x_1, x_2 = i,i+1,i+2
Tnew[x_1] = T(Tbefor[x_0],Tbefor[x_1],Tbefor[x_2],kp,dx,dt)
ti.append(dt*j)
cen.append(Tbefor.tolist())
if j in fre:
plt.plot(x,Tnew,label=str(round(j*dt,6))+"s")
print(j)
Tbefor = Tnew
plt.xlabel("length(um)")
plt.ylabel("Normalaization Thermal")
plt.legend()
plt.xlim(0)
plt.ylim(0,1)
このコードを実行すると以下の画像が描写されます。
空間の分解能が足りていないため、中心が1に正規化されていませんが計算上では問題ないです。
また、時間変化による連続的な熱の拡散は次のように表せます。
max_time = 500
time = [dt*i*10**3 for i in range(max_time)]
length = [dx*i for i in range(Nx)]
cen2 = cen[0:max_time]
def reshape(my_list,n):
le = my_list[0:-1:int(len(my_list)/n)]
return le
sur = go.Surface(y=time,x=length,z=cen2)
fig = go.Figure(sur)
fig.update_layout(scene = dict(
xaxis_title='length (µm)',
yaxis_title='time (ms)',
zaxis_title='Z'),
#width=700,
#margin=dict(r=20, b=10, l=10, t=10)
)
fig.write_html("一次元.html")正弦波的に熱を入射したとき
次に上図のようにsin関数的に熱が入射したときを考えます。さきほどの熱拡散のプログラムに一部改良するだけで大丈夫です。
"""
module import
"""
import matplotlib.pyplot as plt
import numpy as np
import plotly.graph_objects as go
"""
変数定義
"""
r0 = 25 #um -beam radius
xmin,xmax = 0,2000 #um - crystal length
t = 1 #sec - max time
Nt = 10**5 #sample point(time)
Nx = 10**2 #sample point(length)
kp = 16.5/0.683/4*10**6 #熱拡散係数(um^2/s)← 熱拡散係数(m^2/s) = 熱伝導率(W/mK)/(比熱 (J/(g K)) 密度(g/m^3))
TC = 16.5*10**-6 #熱伝導率(W/mK)
dx = (xmax-xmin)/Nx #1step time
dt = (t)/Nt #1step length
Q = 1.5 #power
abs = 0.5 #absorption rate
rep_rate = 1000 #Hz -reptation rate
x = np.linspace(xmin,xmax,Nx) #空間の行列
Tnew = np.zeros(Nx) #初期の温度
print("満たさなければいけない条件(2kpdt<dx**2)",2*kp*dt,"<",dx**2)
"""
数式定義
"""
#ガウス分布
def Gaussian(x,center,r0):
return np.exp(-2*((x-center)/(r0))**2)
#熱伝導方程式
def T(Tx_0,Tx_1,Tx_2,kp,dx,dt):
return Tx_1 + kp*dt*(Tx_0 -2*Tx_1 +Tx_2)/(dx**2)
#正弦波の入熱関数
def sin_funnction(input_Thermal,time,rep_rate):
return input_Thermal*(1-1/2*np.sin(2*np.pi*time*rep_rate))
"""
熱の計算
"""
Tstart = Gaussian(x,x[int(Nx/2)],r0)
T_history = []
Tbefor = Tstart
Thermal = []
time = []
Tbefor[0] = 0
Tbefor[-1] = 0
Gauss_therm = Gaussian(x,x[int(Nx/2)],r0)
for j in range(Nt):
#熱拡散の計算
for i in range(Nx-2):
x_0, x_1, x_2 = i,i+1,i+2
Tnew[x_1] = T(Tbefor[x_0],Tbefor[x_1],Tbefor[x_2],kp,dx,dt)
#熱拡散+加熱
Tnew = Tnew + sin_funnction(Gauss_therm,dt*j,rep_rate)
#時間情報と熱情報の記録
time.append(dt*j)
Thermal.append(Tnew.tolist())
if j % int(Nt/10) == 0:
print("progress:"+str(int(10*(j/Nt+0.11)))+"/10")
#次の計算のために更新
Tbefor = Tnew
"""
3Dで描写
"""
max_time = 500
time = [10000*dt*10**3+dt*i*10**3 for i in range(max_time)]
length = [dx*i for i in range(Nx)]
cen2 = Thermal[0+10000:10000+max_time]
def reshape(my_list,n):
le = my_list[0:-1:int(len(my_list)/n)]
return le
sur = go.Surface(y=time,x=length,z=cen2)
fig = go.Figure(sur)
fig.update_layout(scene = dict(
xaxis_title='length (µm)',
yaxis_title='time (ms)',
zaxis_title='Z'),
#width=700,
#margin=dict(r=20, b=10, l=10, t=10)
)
参考
参考文献:
[1]平田 哲夫,田中 誠,羽田 喜昭,”例題でわかる伝熱工学 第2版”,森北出版(2014)
[2]W. Koechener “Solid-State Laser Engineering”,Springer(2005)
[3]山田啓司 et.al.,” CO2レーザ加工における照射部温度と吸収率”精密工学会誌,65,1,126-130(1999)
[4]J. A. Spenceet.al.,” A review of band structure and material properties of transparent conducting and semiconducting oxides: Ga2O3 ,Al2O3,In2O3,ZnO, SnO2 , CdO, NiO,CuO, and Sc2O3”, Appl Phys Review,9(2022)
参考にしたサイト
【Python】流体シミュレーション : 拡散方程式を実装する
コメント