第5章: テンソル解析と材料科学への応用

import numpy as np # 様々な階数のテンソル scalar = 5.0 # 0階 vector = np.array([1, 2, 3]) # 1階 matrix = np.array([[1, 2], [3, 4]]) # 2階 tensor_3d = np.random.rand(2, 3, 4) # 3階 print("テンソルの例:") print(f"スカラー (0階): {scalar}") print(f" 形状: {np.array(scalar).shape}\n") print(f"ベクトル (1階): {vector}") print(f" 形状: {vector.shape}\n") print(f"行列 (2階):\n{matrix}") print(f" 形状: {matrix.shape}\n") print(f"3階テンソル形状: {tensor_3d.shape}") print(f" 要素数: {tensor_3d.size}\n") # テンソルの転置（軸の入れ替え） tensor_transposed = np.transpose(tensor_3d, (2, 0, 1)) # (2,3,4) -> (4,2,3) print(f"転置後の形状: {tensor_transposed.shape}")

# ベクトルのテンソル積 a = np.array([1, 2, 3]) b = np.array([4, 5]) # 外積（テンソル積） tensor_product = np.outer(a, b) # または np.tensordot(a, b, 0) print("テンソル積:") print(f"a = {a} (形状: {a.shape})") print(f"b = {b} (形状: {b.shape})") print(f"a ⊗ b =\n{tensor_product}") print(f"形状: {tensor_product.shape}\n") # 縮約（トレース） A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) trace = np.trace(A) # 対角成分の和（縮約の例） print(f"行列 A の縮約（トレース）:") print(f"A =\n{A}") print(f"tr(A) = {trace}\n") # Einsteinの縮約記法（np.einsum） # 'ij,jk->ik' は行列積 A @ B A_test = np.array([[1, 2], [3, 4]]) B_test = np.array([[5, 6], [7, 8]]) result_matmul = np.einsum('ij,jk->ik', A_test, B_test) result_direct = A_test @ B_test print("Einsteinの縮約記法:") print(f"A @ B =\n{result_matmul}") print(f"一致? {np.allclose(result_matmul, result_direct)}")

# 任意のテンソル T = np.array([[3, 7, 2], [5, 1, 9], [4, 6, 8]]) # 対称部 T_sym = (T + T.T) / 2 # 反対称部 T_antisym = (T - T.T) / 2 print("対称・反対称分解:") print(f"元のテンソル T:\n{T}\n") print(f"対称部 T_sym = (T + T^T)/2:\n{T_sym}\n") print(f"反対称部 T_antisym = (T - T^T)/2:\n{T_antisym}\n") # 検証: T = T_sym + T_antisym reconstructed = T_sym + T_antisym print(f"再構成 T_sym + T_antisym:\n{reconstructed}") print(f"元と一致? {np.allclose(T, reconstructed)}\n") # 対称性の確認 print(f"T_sym は対称? {np.allclose(T_sym, T_sym.T)}") print(f"T_antisym は反対称? {np.allclose(T_antisym, -T_antisym.T)}")

# 応力テンソル（単位: MPa） stress_tensor = np.array([[100, 30, 0], [30, 60, 0], [0, 0, 20]]) print("応力テンソル解析:") print(f"応力テンソル σ (MPa):\n{stress_tensor}\n") # 対称性の確認 print(f"対称テンソル? {np.allclose(stress_tensor, stress_tensor.T)}\n") # 主応力（固有値）と主軸（固有ベクトル） principal_stresses, principal_axes = np.linalg.eigh(stress_tensor) # 降順にソート idx = principal_stresses.argsort()[::-1] principal_stresses = principal_stresses[idx] principal_axes = principal_axes[:, idx] print("主応力解析:") for i in range(3): print(f"主応力 σ{i+1} = {principal_stresses[i]:.2f} MPa") print(f"主軸方向: {principal_axes[:, i]}\n") # フォン・ミーゼス相当応力 sigma_vm = np.sqrt(((principal_stresses[0] - principal_stresses[1])**2 + (principal_stresses[1] - principal_stresses[2])**2 + (principal_stresses[2] - principal_stresses[0])**2) / 2) print(f"フォン・ミーゼス相当応力: {sigma_vm:.2f} MPa")

# 変位場を定義（簡略版：均一ひずみ） # u_x = ε_xx * x + ε_xy * y # u_y = ε_yx * x + ε_yy * y # 変位勾配テンソル（du_i/dx_j） displacement_gradient = np.array([[0.002, 0.001], # [du_x/dx, du_x/dy] [0.001, 0.003]]) # [du_y/dx, du_y/dy] # ひずみテンソル（対称部） strain_tensor = (displacement_gradient + displacement_gradient.T) / 2 # 回転テンソル（反対称部） rotation_tensor = (displacement_gradient - displacement_gradient.T) / 2 print("ひずみテンソル:") print(f"変位勾配テンソル ∇u:\n{displacement_gradient}\n") print(f"ひずみテンソル ε = (∇u + ∇u^T)/2:\n{strain_tensor}\n") print(f"回転テンソル ω = (∇u - ∇u^T)/2:\n{rotation_tensor}\n") # 体積ひずみ（トレース） volumetric_strain = np.trace(strain_tensor) print(f"体積ひずみ (tr ε): {volumetric_strain:.5f}") # 偏差ひずみテンソル deviatoric_strain = strain_tensor - volumetric_strain/2 * np.eye(2) print(f"偏差ひずみテンソル:\n{deviatoric_strain}")

# 等方性材料の弾性定数 E = 200e9 # ヤング率 (Pa) nu = 0.3 # ポアソン比 # ラメ定数 lambda_lame = (E * nu) / ((1 + nu) * (1 - 2*nu)) mu = E / (2 * (1 + nu)) # せん断弾性率 print("等方性材料の弾性定数:") print(f"ヤング率 E = {E/1e9:.0f} GPa") print(f"ポアソン比 ν = {nu}") print(f"ラメ定数 λ = {lambda_lame/1e9:.2f} GPa") print(f"せん断弾性率 μ = {mu/1e9:.2f} GPa\n") # Voigt記法による弾性テンソル（6×6行列） # [σ_xx, σ_yy, σ_zz, σ_yz, σ_zx, σ_xy]^T = C [ε_xx, ε_yy, ε_zz, 2ε_yz, 2ε_zx, 2ε_xy]^T C_voigt = np.zeros((6, 6)) # 対角成分 C_voigt[0:3, 0:3] = lambda_lame C_voigt += np.diag([2*mu, 2*mu, 2*mu, mu, mu, mu]) print(f"弾性テンソル（Voigt記法、GPa）:") print(np.round(C_voigt/1e9, 2))

# ひずみベクトル（Voigt記法） strain_voigt = np.array([0.001, 0.0005, 0.0003, 0, 0, 0.0002]) # 応力の計算 stress_voigt = C_voigt @ strain_voigt print("\n応力-ひずみ計算:") print(f"ひずみ ε: {strain_voigt}") print(f"応力 σ (MPa): {stress_voigt/1e6}") # テンソル形式に変換 stress_tensor_calc = np.array([[stress_voigt[0], stress_voigt[5], stress_voigt[4]], [stress_voigt[5], stress_voigt[1], stress_voigt[3]], [stress_voigt[4], stress_voigt[3], stress_voigt[2]]]) print(f"\n応力テンソル σ (MPa):\n{stress_tensor_calc/1e6}")

# 立方晶（例: Si）の弾性定数 C11 = 165.7e9 # Pa C12 = 63.9e9 C44 = 79.6e9 # Voigt記法による弾性テンソル C_cubic = np.array([[C11, C12, C12, 0, 0, 0], [C12, C11, C12, 0, 0, 0], [C12, C12, C11, 0, 0, 0], [0, 0, 0, C44, 0, 0], [0, 0, 0, 0, C44, 0], [0, 0, 0, 0, 0, C44]]) print("立方晶の弾性テンソル（Si、GPa）:") print(np.round(C_cubic/1e9, 1)) print(f"\n独立な弾性定数: 3個 (C11, C12, C44)") # 等方性からのずれ（Zener異方性比） A = 2*C44 / (C11 - C12) print(f"\nZener異方性比 A = 2C44/(C11-C12) = {A:.3f}") print(f"(A=1なら完全等方性)") # 体積弾性率と剛性率 K_bulk = (C11 + 2*C12) / 3 G_shear = (C11 - C12 + 3*C44) / 5 print(f"\n平均的性質:") print(f"体積弾性率 K = {K_bulk/1e9:.1f} GPa") print(f"せん断弾性率 G = {G_shear/1e9:.1f} GPa") # ヤング率とポアソン比（等価等方性） E_equivalent = 9*K_bulk*G_shear / (3*K_bulk + G_shear) nu_equivalent = (3*K_bulk - 2*G_shear) / (2*(3*K_bulk + G_shear)) print(f"等価ヤング率 E = {E_equivalent/1e9:.1f} GPa") print(f"等価ポアソン比 ν = {nu_equivalent:.3f}")