情熱！情熱はすべてを解決する！～AIの力で熱伝導のシミュレーションをしてみる
【イントロ】
世の中の物事は、だいたい情熱があれば何でもできます。
そして、熱の流れは固体の物性でも重要です。
固体中の熱の流れは、格子、電子、そして各種素励起により担われます。
一方で、再現性の難しさも指摘されています。
 |
| ©金田一少年の事件簿外伝 犯人たちの事件簿(引用元) |
【実施目的】
熱ホール効果のシミュレーションができればよいのですが、まずは簡単に熱伝導度が試料形状でどのように変化するか、熱伝導のシミュレーションを書いてみることにしました。
形状依存性としては以下の2次元形状のパターンを調べてみました。
- 正方形
- 楔形
- 4つ葉のクローバー型
【手法】
熱伝導の方程式は教科書などに記載されていても、実際の試料中で温度差をつけたときにどのように熱が流れて温度分布が生じるか、シミュレーションする方法がわかりません。
そこで、AIの力を借りました。
Anthropic社のClaudeを利用し、シミュレーションコードを書いてみることにしました。Claudeの有料版機能を利用してコードを作成しました。
原理的には、再現したい論文自体をClaudeに読み込ませ、論文内のシミュレーションを再現することは可能です。ただし、出力されたコードが妥当なのか検証するための知見が足りていないため、その方法は今回は取りませんでした。
作成したコードの実行と結果の可視化はGoogle colabを利用しました。
【結果】
結果です。まずは正方形型から。
 |
| 図１、正方形形状の温度分布 |
見た目上、温度は下部から上部にかけて徐々に高温から低温に変化しているようにみえます。実際、x＝0.5を通る直線上の温度分布を可視化したのが図２です。
 |
| 図２、正方形形状のy軸方向の温度変化 |
直感的にも、滑らかに温度変化しそうなので違和感はない結果になっています。
次に楔形形状です。漢字の「工」みたいな形状の場合です。時代は工学部。
 |
| 図３、楔形形状の温度分布 |
図2と同じようにx＝0.5の線上の温度分布を見てみます。
 |
| 図４、楔形形状のy軸方向の温度変化 |
図4の温度変化を見ると、高温側が形状下部に寄っているようにみえます。
真ん中のくびれによって、熱の流れが詰まっているようなイメージでしょうか？
最後は4つ葉のクローバー形状です。幸運の証を焼きましょう。
 |
| 図５、４つ葉のクローバー形状の温度分布 |
図の4角から中央に向けて45度方向に熱が流れ無い部位があり、形状が4つ葉のクローバーになっています。
温度分布を見ると、高温側、低温側ともに上下方向は1つの葉のなかでは温度変化が見られますが、中央部では変化が一定になっているように見えます。
温度分布を見ると、高温側、低温側ともに上下方向は1つの葉のなかでは温度変化が見られますが、中央部では変化が一定になっているように見えます。
これまでと同じように、x＝0方向の温度変化を可視化してみます。
 |
| 図６、４つ葉のクローバー形状のy軸方向の温度変化 |
図６のように、なんと中央部付近では温度変化が一定になっていることがわかります。
温度差のない左右方向の２つの葉の中で熱の流れが飽和して、中央部では温度変化が一定になっていると解釈できるでしょうか。
【まとめ】
今回は、Claudeを使って2次元形状の熱伝導(熱拡散)による温度分布の変化を可視化してみました。同じ温度差をつけても、サンプル形状により温度分布が変化することが確認できました。試料形状を工夫することで温度分布を操作することに活用できる結果と考えられます。
現実の物質は3次元方向の温度分布もあるため、厚み方向の拡散も考慮にいれる必要があります。また様々な試料を結合して利用する場合は試料ごとの熱伝導率の違いも考慮する必要があります。とはいえ、AIを使えば、簡単なシミュレーションであればすぐに実装できるのは様々な検証を簡単にしてくれたり、勉強した物理や数学の内容を数値計算で確かめることが容易になり学習を深めることができるようになったのは、良い時代の変化です。
一方で、AIが出力する結果が本当に正しいのか、確かめられるスキルが求められるのも確かです。今回の結果が正しいかは、基礎方程式を正しくコーディングできているか利用者が判断できなければなりません。いくら知の高速道路が整備されても、その道を運転する技能を利用者が持っていなければ行けないことは変わらないでしょう。
AIをうまく使って圧倒的成長を遂げたいものです。
【以下、コード】
合ってると思うんだけれども、理解深めないとだめですね（反省）
とはいえ、AIさんすごい。論文も読み込ませられるので、進化を感じました。
時代の先端に食いついていかんとな。
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
# 文字化け対策の設定
plt.rcParams['axes.unicode_minus'] = False # マイナス記号の文字化け防止
plt.rcParams['font.family'] = 'DejaVu Sans' # 英語フォントを明示的に指定
In [ ]:
def solve_heat_equation_2d(nx=50, ny=50, nt=1000, alpha=0.01,
T_hot=100, T_cold=0, boundary_condition='opposite_sides'):
"""
2次元熱拡散方程式を数値的に解く
Parameters:
nx, ny: x, y方向のグリッド数
nt: 時間ステップ数
alpha: 熱拡散係数
T_hot: 高温側の温度
T_cold: 低温側の温度
boundary_condition: 境界条件の種類
"""
# グリッドの設定
dx = 1.0 / (nx - 1)
dy = 1.0 / (ny - 1)
dt = 0.25 * dx * dy / alpha # 安定性のための時間ステップ
# 初期温度分布（全体を低温で初期化）
T = np.ones((ny, nx)) * T_cold
T_new = T.copy()
# 境界条件の設定
if boundary_condition == 'opposite_sides':
# 左右の辺に温度差をかける
T[:, 0] = T_hot # 左辺：高温
T[:, -1] = T_cold # 右辺：低温
# 上下の辺は断熱境界（勾配=0）
elif boundary_condition == 'top_bottom':
# 上下の辺に温度差をかける
T[0, :] = T_hot # 上辺：高温
T[-1, :] = T_cold # 下辺：低温
# 左右の辺は断熱境界
# 時間発展の計算
for n in range(nt):
T_new[1:-1, 1:-1] = (T[1:-1, 1:-1] +
alpha * dt / dx**2 * (T[2:, 1:-1] - 2*T[1:-1, 1:-1] + T[:-2, 1:-1]) +
alpha * dt / dy**2 * (T[1:-1, 2:] - 2*T[1:-1, 1:-1] + T[1:-1, :-2]))
# 境界条件の再適用
if boundary_condition == 'opposite_sides':
T_new[:, 0] = T_hot
T_new[:, -1] = T_cold
# 断熱境界条件（上下）
T_new[0, 1:-1] = T_new[1, 1:-1]
T_new[-1, 1:-1] = T_new[-2, 1:-1]
elif boundary_condition == 'top_bottom':
T_new[0, :] = T_hot
T_new[-1, :] = T_cold
# 断熱境界条件（左右）
T_new[1:-1, 0] = T_new[1:-1, 1]
T_new[1:-1, -1] = T_new[1:-1, -2]
T = T_new.copy()
return T
In [ ]:
def visualize_temperature_gradient():
"""温度勾配を可視化する"""
# パラメータ設定
nx, ny = 50, 50
T_hot = 100 # 高温側 [℃]
T_cold = 0 # 低温側 [℃]
# 2つの境界条件で計算
T1 = solve_heat_equation_2d(nx, ny, nt=2000, T_hot=T_hot, T_cold=T_cold,
boundary_condition='opposite_sides')
T2 = solve_heat_equation_2d(nx, ny, nt=2000, T_hot=T_hot, T_cold=T_cold,
boundary_condition='top_bottom')
# 座標軸の設定
x = np.linspace(0, 1, nx)
y = np.linspace(0, 1, ny)
X, Y = np.meshgrid(x, y)
# 可視化
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# 左右に温度差をかけた場合
im1 = axes[0, 0].contourf(X, Y, T1, levels=20, cmap='coolwarm')
axes[0, 0].set_title('Temperature Distribution (Left-Right Gradient)')
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('y')
fig.colorbar(im1, ax=axes[0, 0], label='Temperature [°C]')
# 等温線表示
cs1 = axes[0, 1].contour(X, Y, T1, levels=10, colors='black', linewidths=0.8)
axes[0, 1].clabel(cs1, inline=True, fontsize=8)
axes[0, 1].set_title('Isotherms (Left-Right Gradient)')
axes[0, 1].set_xlabel('x')
axes[0, 1].set_ylabel('y')
# 上下に温度差をかけた場合
im2 = axes[1, 0].contourf(X, Y, T2, levels=20, cmap='coolwarm')
axes[1, 0].set_title('Temperature Distribution (Top-Bottom Gradient)')
axes[1, 0].set_xlabel('x')
axes[1, 0].set_ylabel('y')
fig.colorbar(im2, ax=axes[1, 0], label='Temperature [°C]')
# 等温線表示
cs2 = axes[1, 1].contour(X, Y, T2, levels=10, colors='black', linewidths=0.8)
axes[1, 1].clabel(cs2, inline=True, fontsize=8)
axes[1, 1].set_title('Isotherms (Top-Bottom Gradient)')
axes[1, 1].set_xlabel('x')
axes[1, 1].set_ylabel('y')
plt.tight_layout()
plt.show()
return T1, T2
In [ ]:
def plot_temperature_profile(T, direction='horizontal'):
"""温度プロファイルをプロット"""
fig, ax = plt.subplots(figsize=(8, 5))
if direction == 'horizontal':
# 中央の水平線での温度プロファイル
middle_row = T.shape[0] // 2
x = np.linspace(0, 1, T.shape[1])
temp_profile = T[middle_row, :]
ax.plot(x, temp_profile, 'b-', linewidth=2, marker='o', markersize=3)
ax.set_xlabel('x coordinate')
ax.set_title('Temperature Profile along Central Horizontal Line')
else:
# 中央の垂直線での温度プロファイル
middle_col = T.shape[1] // 2
y = np.linspace(0, 1, T.shape[0])
temp_profile = T[:, middle_col]
ax.plot(y, temp_profile, 'r-', linewidth=2, marker='s', markersize=3)
ax.set_xlabel('y coordinate')
ax.set_title('Temperature Profile along Central Vertical Line')
ax.set_ylabel('Temperature [°C]')
ax.grid(True, alpha=0.3)
plt.show()
In [ ]:
def create_animation():
"""時間発展のアニメーション"""
nx, ny = 50, 50
nt = 500
T_hot, T_cold = 100, 0
# 初期設定
dx = 1.0 / (nx - 1)
dy = 1.0 / (ny - 1)
alpha = 0.01
dt = 0.25 * dx * dy / alpha
T = np.ones((ny, nx)) * T_cold
T[:, 0] = T_hot # 左辺：高温
T[:, -1] = T_cold # 右辺：低温
x = np.linspace(0, 1, nx)
y = np.linspace(0, 1, ny)
X, Y = np.meshgrid(x, y)
# アニメーション設定
fig, ax = plt.subplots(figsize=(8, 6))
im = ax.contourf(X, Y, T, levels=20, cmap='coolwarm', vmin=T_cold, vmax=T_hot)
cbar = fig.colorbar(im, ax=ax, label='Temperature [°C]')
ax.set_title('Temperature Distribution Evolution')
ax.set_xlabel('x')
ax.set_ylabel('y')
def animate(frame):
nonlocal T
# 時間発展計算
T_new = T.copy()
T_new[1:-1, 1:-1] = (T[1:-1, 1:-1] +
alpha * dt / dx**2 * (T[2:, 1:-1] - 2*T[1:-1, 1:-1] + T[:-2, 1:-1]) +
alpha * dt / dy**2 * (T[1:-1, 2:] - 2*T[1:-1, 1:-1] + T[1:-1, :-2]))
# 境界条件
T_new[:, 0] = T_hot
T_new[:, -1] = T_cold
T_new[0, 1:-1] = T_new[1, 1:-1]
T_new[-1, 1:-1] = T_new[-2, 1:-1]
T = T_new
# プロット更新
ax.clear()
im = ax.contourf(X, Y, T, levels=20, cmap='coolwarm', vmin=T_cold, vmax=T_hot)
ax.set_title(f'Temperature Distribution Evolution (Step: {frame*5})')
ax.set_xlabel('x')
ax.set_ylabel('y')
return im,
anim = FuncAnimation(fig, animate, frames=100, interval=50, blit=False)
plt.show()
return anim
In [ ]:
if __name__ == "__main__":
print("2D Temperature Gradient Visualization")
print("=" * 40)
# 基本的な可視化
T1, T2 = visualize_temperature_gradient()
# 温度プロファイルの表示
print("\nDisplaying temperature profiles...")
plot_temperature_profile(T1, 'horizontal')
plot_temperature_profile(T2, 'vertical')
# # アニメーション（オプション）
# response = input("\nShow time evolution animation? (y/n): ")
# if response.lower() == 'y':
# anim = create_animation()
# print("Animation is being displayed.")
2D Temperature Gradient Visualization ========================================
Displaying temperature profiles...
ここから形状別の可視化
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
# 文字化け対策の設定
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['font.family'] = 'DejaVu Sans'
In [ ]:
def create_constricted_shape(nx=100, ny=100, constriction_type='hourglass'):
"""
形状のマスクを作成
Parameters:
nx, ny: グリッド数
constriction_type: 形状の種類 ('square', 'hourglass', 'neck', 'channel')
Returns:
mask: True=材料内部, False=材料外部
"""
x = np.linspace(0, 1, nx)
y = np.linspace(0, 1, ny)
X, Y = np.meshgrid(x, y)
mask = np.ones((ny, nx), dtype=bool)
if constriction_type == 'square':
# 正方形（全領域が材料）
pass # maskはすでにTrueで初期化されている
elif constriction_type == 'hourglass':
# 砂時計型（上下が広く、中央がくびれ）
center_y = 0.5
center_x = 0.5
for j in range(ny):
y_pos = y[j]
# 中央に近いほど幅が狭くなる
width_factor = 0.2 + 0.6 * (2 * abs(y_pos - center_y))**2
left_bound = center_x - width_factor/2
right_bound = center_x + width_factor/2
for i in range(nx):
x_pos = x[i]
if x_pos < left_bound or x_pos > right_bound:
mask[j, i] = False
elif constriction_type == 'neck':
# ネック型（中央部分だけくびれ）
center_y = 0.5
center_x = 0.5
neck_width = 0.3
neck_height = 0.3
for j in range(ny):
for i in range(nx):
y_pos = y[j]
x_pos = x[i]
# ネック領域（中央部分）
if abs(y_pos - center_y) < neck_height/2:
if abs(x_pos - center_x) > neck_width/2:
mask[j, i] = False
elif constriction_type == 'channel':
# チャンネル型（中央に細い通路）
center_x = 0.5
channel_width = 0.2
for j in range(ny):
for i in range(nx):
x_pos = x[i]
y_pos = y[j]
# 上下の領域は広い
if y_pos < 0.3 or y_pos > 0.7:
continue # 全幅使用
else:
# 中央部分は狭いチャンネル
if abs(x_pos - center_x) > channel_width/2:
mask[j, i] = False
return mask
In [ ]:
def solve_heat_equation_constricted(nx=100, ny=100, nt=2000, alpha=0.01,
T_hot=100, T_cold=0, constriction_type='hourglass'):
"""
様々な形状での熱拡散方程式を解く
"""
# 形状マスクの作成
mask = create_constricted_shape(nx, ny, constriction_type)
# グリッドの設定
dx = 1.0 / (nx - 1)
dy = 1.0 / (ny - 1)
dt = 0.1 * dx * dy / alpha # 安定性のための時間ステップ
# 初期温度分布
T = np.ones((ny, nx)) * T_cold
T_new = T.copy()
# 境界条件：上辺を高温、下辺を低温に設定
for i in range(nx):
if mask[0, i]: # 上辺（y=0, 配列の最初の行）
T[0, i] = T_hot
if mask[-1, i]: # 下辺（y=1, 配列の最後の行）
T[-1, i] = T_cold
# 材料外部は温度計算しない（NaNで表示）
T[~mask] = np.nan
T_new[~mask] = np.nan
# 時間発展計算
for n in range(nt):
for j in range(1, ny-1):
for i in range(1, nx-1):
if mask[j, i]: # 材料内部のみ計算
# 隣接セルが材料内部かチェック
neighbors = []
if mask[j+1, i]: # 下
neighbors.append(T[j+1, i])
else:
neighbors.append(T[j, i]) # 断熱境界
if mask[j-1, i]: # 上
neighbors.append(T[j-1, i])
else:
neighbors.append(T[j, i])
if mask[j, i+1]: # 右
neighbors.append(T[j, i+1])
else:
neighbors.append(T[j, i])
if mask[j, i-1]: # 左
neighbors.append(T[j, i-1])
else:
neighbors.append(T[j, i])
# 熱拡散方程式
T_new[j, i] = (T[j, i] +
alpha * dt / dx**2 * (neighbors[3] - 2*T[j, i] + neighbors[2]) +
alpha * dt / dy**2 * (neighbors[0] - 2*T[j, i] + neighbors[1]))
# 境界条件の再適用
for i in range(nx):
if mask[0, i]: # 上辺
T_new[0, i] = T_hot
if mask[-1, i]: # 下辺
T_new[-1, i] = T_cold
T = T_new.copy()
# 収束チェック（オプション）
if n % 500 == 0:
print(f"Step {n}/{nt}")
return T, mask
In [ ]:
def visualize_single_shape(shape_type='neck'):
"""指定した形状での温度勾配を可視化"""
# パラメータ設定
nx, ny = 100, 100
T_hot = 100
T_cold = 0
# 形状名の辞書
shape_names = {
'square': 'Square Shape',
'hourglass': 'Hourglass Shape',
'neck': 'Neck Shape',
'channel': 'Channel Shape'
}
if shape_type not in shape_names:
print(f"Invalid shape type: {shape_type}")
print(f"Available shapes: {list(shape_names.keys())}")
return None
# 指定された形状で計算
T, mask = solve_heat_equation_constricted(nx, ny, nt=3000,
T_hot=T_hot, T_cold=T_cold,
constriction_type=shape_type)
x = np.linspace(0, 1, nx)
y = np.linspace(0, 1, ny)
X, Y = np.meshgrid(x, y)
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 温度分布
im1 = axes[0].contourf(X, Y, T, levels=20, cmap='coolwarm',
vmin=T_cold, vmax=T_hot)
axes[0].set_title(f'Temperature Distribution\n{shape_names[shape_type]}', fontsize=14)
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_aspect('equal')
# カラーバーを追加
plt.colorbar(im1, ax=axes[0], label='Temperature [°C]')
# 等温線と形状境界
T_masked = T.copy()
T_masked[~mask] = np.nan
axes[1].contour(X, Y, T_masked, levels=15, colors='black', linewidths=0.8)
axes[1].contourf(X, Y, mask.astype(float), levels=[0, 0.5, 1],
colors=['white', 'lightgray'], alpha=0.3)
axes[1].set_title(f'Isotherms and Shape\n{shape_names[shape_type]}', fontsize=14)
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].set_aspect('equal')
plt.tight_layout()
plt.show()
return T, mask
In [ ]:
def visualize_all_shapes():
"""全ての形状での温度勾配を可視化"""
# パラメータ設定
nx, ny = 100, 100
T_hot = 100
T_cold = 0
# 4つの形状タイプで計算
shapes = ['square', 'hourglass', 'neck', 'channel']
shape_names = ['Square Shape', 'Hourglass Shape', 'Neck Shape', 'Channel Shape']
fig, axes = plt.subplots(2, 4, figsize=(20, 10))
for idx, (shape_type, shape_name) in enumerate(zip(shapes, shape_names)):
T, mask = solve_heat_equation_constricted(nx, ny, nt=3000,
T_hot=T_hot, T_cold=T_cold,
constriction_type=shape_type)
x = np.linspace(0, 1, nx)
y = np.linspace(0, 1, ny)
X, Y = np.meshgrid(x, y)
# 温度分布
im1 = axes[0, idx].contourf(X, Y, T, levels=20, cmap='coolwarm',
vmin=T_cold, vmax=T_hot)
axes[0, idx].set_title(f'Temperature Distribution\n{shape_name}')
axes[0, idx].set_xlabel('x')
axes[0, idx].set_ylabel('y')
axes[0, idx].set_aspect('equal')
# カラーバーを追加
plt.colorbar(im1, ax=axes[0, idx], label='Temperature [°C]', shrink=0.8)
# 等温線と形状境界
# 材料領域のみを表示するため、マスクされた領域を白で塗りつぶし
T_masked = T.copy()
T_masked[~mask] = np.nan
axes[1, idx].contour(X, Y, T_masked, levels=15, colors='black', linewidths=0.8)
axes[1, idx].contourf(X, Y, mask.astype(float), levels=[0, 0.5, 1],
colors=['white', 'lightgray'], alpha=0.3)
axes[1, idx].set_title(f'Isotherms and Shape\n{shape_name}')
axes[1, idx].set_xlabel('x')
axes[1, idx].set_ylabel('y')
axes[1, idx].set_aspect('equal')
plt.tight_layout()
plt.show()
return shapes
In [ ]:
def plot_temperature_along_centerline(shape_type='hourglass'):
"""中心線に沿った温度プロファイル"""
nx, ny = 100, 100
T, mask = solve_heat_equation_constricted(nx, ny, nt=3000, constriction_type=shape_type)
# 中心線（x=0.5）での温度プロファイル
center_col = nx // 2
y = np.linspace(0, 1, ny)
# 材料内部のみのデータを抽出
valid_indices = mask[:, center_col]
y_valid = y[valid_indices]
T_valid = T[valid_indices, center_col]
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(y_valid, T_valid, 'ro-', linewidth=2, markersize=4)
ax.set_xlabel('y coordinate')
ax.set_ylabel('Temperature [°C]')
ax.set_title(f'Temperature Profile along Centerline ({shape_type.capitalize()} Shape)')
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 1)
plt.show()
In [ ]:
if __name__ == "__main__":
print("Shape Temperature Gradient Visualization")
print("=" * 40)
# ネック型の可視化
print("Visualizing Neck Shape...")
T_neck, mask_neck = visualize_single_shape('neck')
# 中心線温度プロファイル（ネック型）
print("\nDisplaying centerline temperature profile for Neck shape...")
plot_temperature_along_centerline('neck')
# # 他の形状も表示したい場合のオプション
# print("\nOther available shapes: square, hourglass, channel")
# print("Use visualize_single_shape('shape_name') to visualize specific shape")
# print("Use visualize_all_shapes() to compare all shapes")
Shape Temperature Gradient Visualization ======================================== Visualizing Neck Shape... Step 0/3000 Step 500/3000 Step 1000/3000 Step 1500/3000 Step 2000/3000 Step 2500/3000
Displaying centerline temperature profile for Neck shape... Step 0/3000 Step 500/3000 Step 1000/3000 Step 1500/3000 Step 2000/3000 Step 2500/3000
In [ ]:
# 中心線温度プロファイル（ネック型）
print("\nDisplaying centerline temperature profile for Neck shape...")
plot_temperature_along_centerline('neck')
Displaying centerline temperature profile for Neck shape... Step 0/3000 Step 500/3000 Step 1000/3000 Step 1500/3000 Step 2000/3000 Step 2500/3000
In [ ]:
4つばのクローバー形状
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
# 文字化け対策の設定
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['font.family'] = 'DejaVu Sans'
In [ ]:
def create_cross_shape(nx=100, ny=100, rho=0.3):
"""
十字型形状のマスクを作成
大きな正方形から対角線上（中央正方形以外）を削除した形状
Parameters:
nx, ny: グリッド数
rho: 中央正方形のサイズ比（中央正方形の一辺/全体の一辺）
Returns:
mask: True=材料内部, False=削り取られた部分（空洞）
"""
x = np.linspace(0, 1, nx)
y = np.linspace(0, 1, ny)
# 初期化：全体をTrue（大きな正方形A）
mask = np.ones((ny, nx), dtype=bool)
# 中央正方形Bの範囲
center = 0.5
half_b = rho / 2
for j in range(ny):
for i in range(nx):
x_pos = x[i]
y_pos = y[j]
# 中央正方形B内部かチェック
in_center_square = (center - half_b <= x_pos <= center + half_b and
center - half_b <= y_pos <= center + half_b)
# 中央正方形内部なら必ず材料として保持
if in_center_square:
mask[j, i] = True
else:
# 中央正方形外部で対角線上かチェック
# 対角線の幅（線の太さ）
line_width = 1.0 / nx # グリッド1つ分の幅
# 対角線1: y = x
on_diag1 = abs(y_pos - x_pos) <= line_width
# 対角線2: y = 1 - x
on_diag2 = abs(y_pos - (1 - x_pos)) <= line_width
# 対角線上で、かつ中央正方形外部なら削除
if on_diag1 or on_diag2:
mask[j, i] = False
else:
# 対角線上でなければ材料として保持
mask[j, i] = True
return mask
In [ ]:
def solve_heat_equation_cross(nx=100, ny=100, nt=3000, alpha=0.01,
T_hot=100, T_cold=0, rho=0.3):
"""
十字型形状での熱拡散方程式を解く
"""
# 形状マスクの作成
mask = create_cross_shape(nx, ny, rho)
# グリッドの設定
dx = 1.0 / (nx - 1)
dy = 1.0 / (ny - 1)
dt = 0.05 * dx * dy / alpha # 安定性のための時間ステップ
# 初期温度分布
T = np.ones((ny, nx)) * ((T_hot + T_cold) / 2) # 中間温度で初期化
T_new = T.copy()
# 境界条件設定
center = 0.5
half_d = rho / 2
for j in range(ny):
for i in range(nx):
if mask[j, i]: # 材料内部のみ
x_pos = i / (nx - 1)
y_pos = j / (ny - 1)
# 上辺境界（heater）
if j == 0:
T[j, i] = T_hot
# 下辺境界（heat sink）
elif j == ny - 1:
T[j, i] = T_cold
# 左辺境界（thermometer - 断熱）
elif i == 0:
pass # 断熱境界（後で隣接セルと同じ温度に設定）
# 右辺境界（thermometer - 断熱）
elif i == nx - 1:
pass # 断熱境界
# 材料外部は計算しない
T[~mask] = np.nan
T_new[~mask] = np.nan
# 時間発展計算
for n in range(nt):
for j in range(1, ny-1):
for i in range(1, nx-1):
if mask[j, i]: # 材料内部のみ計算
# 隣接セルの温度を取得（材料外部の場合は現在の温度を使用）
T_up = T[j-1, i] if mask[j-1, i] else T[j, i]
T_down = T[j+1, i] if mask[j+1, i] else T[j, i]
T_left = T[j, i-1] if mask[j, i-1] else T[j, i]
T_right = T[j, i+1] if mask[j, i+1] else T[j, i]
# 熱拡散方程式
T_new[j, i] = (T[j, i] +
alpha * dt / dx**2 * (T_left - 2*T[j, i] + T_right) +
alpha * dt / dy**2 * (T_up - 2*T[j, i] + T_down))
# 境界条件の再適用
for j in range(ny):
for i in range(nx):
if mask[j, i]:
# 上辺（heater）
if j == 0:
T_new[j, i] = T_hot
# 下辺（heat sink）
elif j == ny - 1:
T_new[j, i] = T_cold
# 左右辺（断熱境界）
elif i == 0 and j > 0 and j < ny-1:
if mask[j, 1]:
T_new[j, i] = T_new[j, 1]
elif i == nx-1 and j > 0 and j < ny-1:
if mask[j, nx-2]:
T_new[j, i] = T_new[j, nx-2]
T = T_new.copy()
# 収束チェック
if n % 1000 == 0:
print(f"Step {n}/{nt}")
return T, mask
In [ ]:
def visualize_cross_shape(rho=0.3):
"""十字型形状での温度勾配を可視化"""
# パラメータ設定
nx, ny = 120, 120
T_hot = 100
T_cold = 0
print(f"Computing temperature distribution for cross shape (ρ = {rho:.2f})...")
T, mask = solve_heat_equation_cross(nx, ny, nt=4000,
T_hot=T_hot, T_cold=T_cold, rho=rho)
x = np.linspace(0, 1, nx)
y = np.linspace(0, 1, ny)
X, Y = np.meshgrid(x, y)
# 可視化
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
# 1. 形状の表示
shape_display = mask.astype(float)
shape_display[~mask] = 0 # 材料外部は0
im0 = axes[0].imshow(shape_display, cmap='gray', origin='lower', extent=[0, 1, 0, 1])
axes[0].set_title(f'Cross Shape Structure\n(ρ = d/D = {rho:.2f})', fontsize=14)
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].text(0.5, 0.9, 'HEATER', ha='center', va='center',
color='red', fontsize=12, weight='bold')
axes[0].text(0.5, 0.1, 'HEAT SINK', ha='center', va='center',
color='blue', fontsize=12, weight='bold')
axes[0].text(0.05, 0.5, 'THERMO', ha='center', va='center',
rotation=90, color='green', fontsize=10, weight='bold')
axes[0].text(0.95, 0.5, 'THERMO', ha='center', va='center',
rotation=90, color='green', fontsize=10, weight='bold')
# 2. 温度分布
im1 = axes[1].contourf(X, Y, T, levels=20, cmap='coolwarm',
vmin=T_cold, vmax=T_hot)
axes[1].set_title(f'Temperature Distribution\nCross Shape (ρ = {rho:.2f})', fontsize=14)
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
axes[1].set_aspect('equal')
plt.colorbar(im1, ax=axes[1], label='Temperature [°C]')
# 3. 等温線
T_masked = T.copy()
T_masked[~mask] = np.nan
axes[2].contour(X, Y, T_masked, levels=15, colors='black', linewidths=0.8)
axes[2].contourf(X, Y, mask.astype(float), levels=[0, 0.5, 1],
colors=['white', 'lightgray'], alpha=0.3)
axes[2].set_title(f'Isotherms\nCross Shape (ρ = {rho:.2f})', fontsize=14)
axes[2].set_xlabel('x')
axes[2].set_ylabel('y')
axes[2].set_aspect('equal')
plt.tight_layout()
plt.show()
return T, mask
In [ ]:
def plot_temperature_profiles_cross(rho=0.3):
"""十字型形状での温度プロファイル"""
nx, ny = 120, 120
T, mask = solve_heat_equation_cross(nx, ny, nt=4000, rho=rho)
fig, axes = plt.subplots(1, 2, figsize=(15, 6))
# 垂直方向プロファイル（中央）
center_col = nx // 2
y = np.linspace(0, 1, ny)
valid_indices = mask[:, center_col]
y_valid = y[valid_indices]
T_valid = T[valid_indices, center_col]
axes[0].plot(y_valid, T_valid, 'ro-', linewidth=2, markersize=4)
axes[0].set_xlabel('y coordinate')
axes[0].set_ylabel('Temperature [°C]')
axes[0].set_title(f'Vertical Temperature Profile\n(ρ = {rho:.2f})')
axes[0].grid(True, alpha=0.3)
axes[0].set_xlim(0, 1)
# 水平方向プロファイル（中央）
center_row = ny // 2
x = np.linspace(0, 1, nx)
valid_indices = mask[center_row, :]
x_valid = x[valid_indices]
T_valid = T[center_row, valid_indices]
axes[1].plot(x_valid, T_valid, 'bo-', linewidth=2, markersize=4)
axes[1].set_xlabel('x coordinate')
axes[1].set_ylabel('Temperature [°C]')
axes[1].set_title(f'Horizontal Temperature Profile\n(ρ = {rho:.2f})')
axes[1].grid(True, alpha=0.3)
axes[1].set_xlim(0, 1)
plt.tight_layout()
plt.show()
In [ ]:
def compare_rho_values():
"""異なるρ値での比較"""
rho_values = [0.2, 0.3, 0.4, 0.5]
fig, axes = plt.subplots(2, 4, figsize=(20, 10))
for idx, rho in enumerate(rho_values):
print(f"Computing for ρ = {rho}")
nx, ny = 100, 100
T, mask = solve_heat_equation_cross(nx, ny, nt=3000, rho=rho)
x = np.linspace(0, 1, nx)
y = np.linspace(0, 1, ny)
X, Y = np.meshgrid(x, y)
# 温度分布
im1 = axes[0, idx].contourf(X, Y, T, levels=20, cmap='coolwarm',
vmin=0, vmax=100)
axes[0, idx].set_title(f'Temperature (ρ = {rho})')
axes[0, idx].set_xlabel('x')
axes[0, idx].set_ylabel('y')
axes[0, idx].set_aspect('equal')
# 等温線
T_masked = T.copy()
T_masked[~mask] = np.nan
axes[1, idx].contour(X, Y, T_masked, levels=10, colors='black', linewidths=0.8)
axes[1, idx].contourf(X, Y, mask.astype(float), levels=[0, 0.5, 1],
colors=['white', 'lightgray'], alpha=0.3)
axes[1, idx].set_title(f'Isotherms (ρ = {rho})')
axes[1, idx].set_xlabel('x')
axes[1, idx].set_ylabel('y')
axes[1, idx].set_aspect('equal')
plt.tight_layout()
plt.show()
In [ ]:
if __name__ == "__main__":
print("Cross-Shaped Sample Temperature Gradient Visualization")
print("=" * 55)
# デフォルトのρ値で可視化
rho_default = 0.3
print(f"\nVisualizing cross shape with ρ = {rho_default}...")
T, mask = visualize_cross_shape(rho_default)
# 温度プロファイル表示
print(f"\nDisplaying temperature profiles...")
plot_temperature_profiles_cross(rho_default)
# # 異なるρ値での比較
# print(f"\nComparing different ρ values...")
# compare_rho_values()
# print("\nAvailable functions:")
# print("- visualize_cross_shape(rho): Single ρ value visualization")
# print("- plot_temperature_profiles_cross(rho): Temperature profiles")
# print("- compare_rho_values(): Compare multiple ρ values")
Cross-Shaped Sample Temperature Gradient Visualization ======================================================= Visualizing cross shape with ρ = 0.3... Computing temperature distribution for cross shape (ρ = 0.30)... Step 0/4000 Step 1000/4000 Step 2000/4000 Step 3000/4000
Displaying temperature profiles... Step 0/4000 Step 1000/4000 Step 2000/4000 Step 3000/4000
In [ ]:
コメント
コメントを投稿