import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import g
# Viscosity
def viscosity_plot():
mu = 0.001 # Viscosity (Pa·s), e.g., viscosity of water
u = np.linspace(0, 1, 100)
shear_stress = mu * u
plt.figure(figsize=(10, 6))
plt.plot(u, shear_stress, label=f'Viscosity μ = {mu} Pa·s')
plt.xlabel('Shear Rate (s⁻¹)')
plt.ylabel('Shear Stress (Pa)')
plt.title('Viscosity Plot')
plt.legend()
plt.grid(True)
plt.show()
# Density
def density_plot():
densities = {'Water': 1000, 'Oil': 800, 'Air': 1.225} # kg/m³
substances = list(densities.keys())
values = list(densities.values())
plt.figure(figsize=(10, 6))
plt.bar(substances, values, color=['blue', 'orange', 'green'])
plt.xlabel('Substance')
plt.ylabel('Density (kg/m³)')
plt.title('Density of Various Substances')
plt.grid(axis='y')
plt.show()
# Pressure
def pressure_plot():
h = np.linspace(0, 10, 100) # Height (m)
rho = 1000 # Density (kg/m³)
g = 9.81 # Gravitational Acceleration (m/s²)
pressure = rho * g * h
plt.figure(figsize=(10, 6))
plt.plot(h, pressure, label='Pressure')
plt.xlabel('Height (m)')
plt.ylabel('Pressure (Pa)')
plt.title('Pressure Plot')
plt.legend()
plt.grid(True)
plt.show()
# Surface Tension and Capillary Action
def capillary_action_plot():
surface_tension = 0.0728 # Surface Tension (N/m), e.g., water
radius = np.linspace(0.1, 1, 100)
rho = 1000 # Density of water (kg/m³)
height = (2 * surface_tension) / (rho * g * radius) # Capillary Rise Height
plt.figure(figsize=(10, 6))
plt.plot(radius, height, label='Capillary Rise Height')
plt.xlabel('Tube Radius (m)')
plt.ylabel('Height (m)')
plt.title('Capillary Action Plot')
plt.legend()
plt.grid(True)
plt.show()
# Display plots
viscosity_plot()
density_plot()
pressure_plot()
capillary_action_plot()
import numpy as np
import matplotlib.pyplot as plt
# Constants
density_water = 1000 # Density of water [kg/m^3]
pressure_range = np.linspace(0, 1e6, 500) # Pressure range [Pa]
# Calculate density (example: density depends on pressure)
density_fluid = 900 + 0.1 * pressure_range # [kg/m^3]
# Calculate specific gravity
specific_gravity = density_fluid / density_water
# Calculate compressibility (beta is derived from density change rate)
compressibility = (density_fluid - 900) / (pressure_range + 1) # [1/Pa]
# Plot density and specific gravity
fig, ax1 = plt.subplots()
ax1.set_xlabel('Pressure [Pa]')
ax1.set_ylabel('Density [kg/m^3]', color='tab:blue')
ax1.plot(pressure_range, density_fluid, color='tab:blue', label='Density')
ax1.tick_params(axis='y', labelcolor='tab:blue')
ax2 = ax1.twinx()
ax2.set_ylabel('Specific Gravity', color='tab:orange')
ax2.plot(pressure_range, specific_gravity, color='tab:orange', label='Specific Gravity')
ax2.tick_params(axis='y', labelcolor='tab:orange')
fig.tight_layout()
plt.title('Density and Specific Gravity vs Pressure')
plt.show()
# Plot compressibility
fig, ax = plt.subplots()
ax.plot(pressure_range, compressibility, label='Compressibility')
ax.set_xlabel('Pressure [Pa]')
ax.set_ylabel('Compressibility [1/Pa]')
ax.set_title('Compressibility vs Pressure')
ax.legend()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Constants
g = 9.81 # Acceleration due to gravity [m/s^2]
density_water = 1000 # Density of water [kg/m^3]
density_air = 1.225 # Density of air [kg/m^3]
height_range = np.linspace(0, 10, 500) # Water column height range [m]
# 2.1 Hydrostatic Pressure Calculation
def hydrostatic_pressure(height, density=density_water):
return density * g * height
# 2.2 Buoyancy Calculation
def buoyant_force(volume, density=density_water):
return density * g * volume
# 2.3 Pressure Difference Calculation
def pressure_difference(height):
return hydrostatic_pressure(height) - density_air * g * height
# Hydrostatic pressure and buoyancy calculations
pressure_water = hydrostatic_pressure(height_range)
buoyant_force_water = buoyant_force(height_range * 1) # Assume 1 m^3 object for buoyancy calculation
pressure_diff = pressure_difference(height_range)
# Plot Hydrostatic Pressure
fig, ax1 = plt.subplots()
ax1.set_xlabel('Height of Water Column [m]')
ax1.set_ylabel('Hydrostatic Pressure [Pa]', color='tab:blue')
ax1.plot(height_range, pressure_water, color='tab:blue', label='Hydrostatic Pressure')
ax1.tick_params(axis='y', labelcolor='tab:blue')
# Plot Buoyant Force
ax2 = ax1.twinx()
ax2.set_ylabel('Buoyant Force [N]', color='tab:orange')
ax2.plot(height_range, buoyant_force_water, color='tab:orange', label='Buoyant Force')
ax2.tick_params(axis='y', labelcolor='tab:orange')
fig.tight_layout()
plt.title('Hydrostatic Pressure and Buoyant Force vs Height')
plt.show()
# Plot Pressure Difference
fig, ax = plt.subplots()
ax.plot(height_range, pressure_diff, label='Pressure Difference')
ax.set_xlabel('Height of Water Column [m]')
ax.set_ylabel('Pressure Difference [Pa]')
ax.set_title('Pressure Difference between Water and Air vs Height')
ax.legend()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Constants
area1 = 1 # Cross-sectional area at point 1 [m^2]
velocity1 = 2 # Velocity at point 1 [m/s]
# Define cross-sectional area range for point 2
area2_range = np.linspace(0.1, 2, 500) # [m^2]
# Calculate velocity at point 2 using continuity equation
velocity2_range = (area1 * velocity1) / area2_range
# Plot velocity vs. cross-sectional area
plt.figure(figsize=(8, 6))
plt.plot(area2_range, velocity2_range, label='Velocity at point 2')
plt.xlabel('Cross-sectional Area at point 2 [m^2]')
plt.ylabel('Velocity at point 2 [m/s]')
plt.title('Velocity vs Cross-sectional Area using Continuity Equation')
plt.legend()
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Parameters for laminar flow in a pipe
radius = 1 # Radius of the pipe [m]
viscosity = 0.001 # Viscosity [Pa.s]
pressure_drop = 1000 # Pressure drop over the length [Pa]
# Define radial position from center to edge
r = np.linspace(0, radius, 500)
# Calculate velocity profile using simplified Navier-Stokes for laminar flow
velocity_profile = (pressure_drop * (radius**2 - r**2)) / (4 * viscosity * radius)
# Plot velocity profile
plt.figure(figsize=(8, 6))
plt.plot(r, velocity_profile, label='Velocity Profile')
plt.xlabel('Radial Position [m]')
plt.ylabel('Velocity [m/s]')
plt.title('Velocity Profile for Laminar Flow in a Pipe')
plt.legend()
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# ドメインサイズとグリッドの設定
nx, ny = 50, 50
dx, dy = 1.0, 1.0
nt = 1000
rho = 1.0
mu = 0.1
# 初期化
u = np.zeros((nx, ny))
v = np.zeros((nx, ny))
p = np.zeros((nx, ny))
# ナビエ-ストークス方程式の簡略化された解法
def solve_navier_stokes(u, v, p, mu, dx, dy, nt):
for _ in range(nt):
u_old = u.copy()
v_old = v.copy()
# 圧力ポアソン方程式を解く（簡略化された例）
p[1:-1, 1:-1] = (1/4) * (p[2:, 1:-1] + p[:-2, 1:-1] + p[1:-1, 2:] + p[1:-1, :-2])
# x方向の速度更新
u[1:-1, 1:-1] = u_old[1:-1, 1:-1] + (mu * (
(u_old[2:, 1:-1] - 2 * u_old[1:-1, 1:-1] + u_old[:-2, 1:-1]) / dx**2 +
(u_old[1:-1, 2:] - 2 * u_old[1:-1, 1:-1] + u_old[1:-1, :-2]) / dy**2
))
# y方向の速度更新
v[1:-1, 1:-1] = v_old[1:-1, 1:-1] + (mu * (
(v_old[2:, 1:-1] - 2 * v_old[1:-1, 1:-1] + v_old[:-2, 1:-1]) / dx**2 +
(v_old[1:-1, 2:] - 2 * v_old[1:-1, 1:-1] + v_old[1:-1, :-2]) / dy**2
))
return u, v, p
# 解の計算
u, v, p = solve_navier_stokes(u, v, p, mu, dx, dy, nt)
# 結果のプロット
plt.figure(figsize=(12, 6))
plt.subplot(1, 3, 1)
plt.title('Velocity field (u)')
plt.imshow(u.T, origin='lower', cmap='jet')
plt.colorbar()
plt.subplot(1, 3, 2)
plt.title('Velocity field (v)')
plt.imshow(v.T, origin='lower', cmap='jet')
plt.colorbar()
plt.subplot(1, 3, 3)
plt.title('Pressure field')
plt.imshow(p.T, origin='lower', cmap='jet')
plt.colorbar()
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# パラメータ設定
rho = 1.0 # 流体の密度 (kg/m^3)
g = 9.81 # 重力加速度 (m/s^2)
h1 = 10.0 # 位置1の高さ (m)
h2 = 0.0 # 位置2の高さ (m)
p1 = 100000 # 位置1の圧力 (Pa)
# 高さと速度の範囲を設定
heights = np.linspace(0, 10, 100)
velocities = np.zeros_like(heights)
pressures = np.zeros_like(heights)
# 高さに基づく速度と圧力の計算
for i, h in enumerate(heights):
p2 = p1 + rho * g * (h1 - h)
velocities[i] = np.sqrt(2 * (p1 - p2) / rho) # 速度の計算
pressures[i] = p2
# 結果のプロット
plt.figure(figsize=(12, 6))
plt.subplot(1, 3, 1)
plt.title('Velocity vs Height')
plt.plot(heights, velocities, label='Velocity')
plt.xlabel('Height (m)')
plt.ylabel('Velocity (m/s)')
plt.legend()
plt.subplot(1, 3, 2)
plt.title('Pressure vs Height')
plt.plot(heights, pressures, label='Pressure', color='r')
plt.xlabel('Height (m)')
plt.ylabel('Pressure (Pa)')
plt.legend()
plt.subplot(1, 3, 3)
plt.title('Total Energy vs Height')
total_energy = 0.5 * rho * velocities**2 + rho * g * heights
plt.plot(heights, total_energy, label='Total Energy', color='g')
plt.xlabel('Height (m)')
plt.ylabel('Total Energy (J/m^3)')
plt.legend()
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# グリッドの設定
x = np.linspace(0, 10, 100)
y = np.linspace(0, 10, 100)
X, Y = np.meshgrid(x, y)
# 速度場の定義（例として2次元のポテンシャル流れを使用）
U = -np.sin(np.pi * X) * np.cos(np.pi * Y)
V = np.cos(np.pi * X) * np.sin(np.pi * Y)
# 流線の描画
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.title('Streamlines')
plt.streamplot(X, Y, U, V, color='b', linewidth=1.5)
plt.xlabel('X')
plt.ylabel('Y')
# ベクトル場の描画
plt.subplot(1, 2, 2)
plt.title('Velocity Field')
plt.quiver(X, Y, U, V, scale=50)
plt.xlabel('X')
plt.ylabel('Y')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# パラメータ設定
L = 1.0 # 代表長さ (m)
u_inf = 1.0 # 流速 (m/s)
nu = 1e-6 # 動粘性係数 (m^2/s)
x = np.linspace(0.01, L, 100) # 流れ方向の範囲
# 境界層の厚さの計算
# デルタ* = 5 * sqrt(nu * x / u_inf) という近似式を使用
delta_star = 5 * np.sqrt(nu * x / u_inf)
# 速度プロファイルの計算
y = np.linspace(0, 2*delta_star.max(), 100)
U = np.zeros((len(x), len(y)))
for i, x_i in enumerate(x):
U[i, :] = u_inf * (1 - np.exp(-y / (5 * np.sqrt(nu * x_i / u_inf))))
# プロット
plt.figure(figsize=(12, 6))
# 境界層の厚さのプロット
plt.subplot(1, 2, 1)
plt.title('Boundary Layer Thickness')
plt.plot(x, delta_star, label='Boundary Layer Thickness')
plt.xlabel('Distance from Leading Edge (m)')
plt.ylabel('Boundary Layer Thickness (m)')
plt.legend()
# 速度プロファイルのプロット
X, Y = np.meshgrid(x, y)
plt.subplot(1, 2, 2)
plt.title('Velocity Profile')
plt.contourf(X, Y, U.T, levels=20, cmap='viridis')
plt.colorbar(label='Velocity (m/s)')
plt.xlabel('Distance from Leading Edge (m)')
plt.ylabel('Height (m)')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
# ポンプの設計パラメータ
def pump_head(Q, Q_max, H_max):
return H_max * (1 - (Q / Q_max)**2)
# ファンの設計パラメータ
def fan_pressure(Q, Q_max, P_max):
return P_max * (1 - (Q / Q_max)**2)
# ポンプの流量と揚程の関係をプロット
Q_pump = np.linspace(0, 100, 100)
H_pump = pump_head(Q_pump, 100, 50)
# ファンの風量と静圧の関係をプロット
Q_fan = np.linspace(0, 100, 100)
P_fan = fan_pressure(Q_fan, 100, 500)
plt.figure(figsize=(12, 6))
# ポンプのプロット
plt.subplot(1, 2, 1)
plt.title('Pump Flow Rate vs. Head')
plt.plot(Q_pump, H_pump, label='Pump Head')
plt.xlabel('Flow Rate (m³/h)')
plt.ylabel('Head (m)')
plt.legend()
# ファンのプロット
plt.subplot(1, 2, 2)
plt.title('Fan Flow Rate vs. Pressure')
plt.plot(Q_fan, P_fan, label='Fan Pressure', color='r')
plt.xlabel('Flow Rate (m³/h)')
plt.ylabel('Pressure (Pa)')
plt.legend()
plt.tight_layout()
plt.show()
# 流れの制御と最適化（簡単な例）
def objective(x):
# ここでは単純に流量と圧力のバランスを最適化
Q_opt, P_opt = x
return (Q_opt - 50)**2 + (P_opt - 250)**2 # 最適化目的関数
# 制約条件
def constraint1(x):
return x[0] - 20 # 最小流量制約
def constraint2(x):
return 80 - x[0] # 最大流量制約
def constraint3(x):
return x[1] - 100 # 最小圧力制約
def constraint4(x):
return 400 - x[1] # 最大圧力制約
x0 = [30, 200] # 初期推定値
constraints = [{'type': 'ineq', 'fun': constraint1},
{'type': 'ineq', 'fun': constraint2},
{'type': 'ineq', 'fun': constraint3},
{'type': 'ineq', 'fun': constraint4}]
result = minimize(objective, x0, constraints=constraints)
print("Optimal flow rate:", result.x[0])
print("Optimal pressure:", result.x[1])
import numpy as np
import matplotlib.pyplot as plt
# メッシュの生成
def generate_mesh(xmin, xmax, ymin, ymax, nx, ny):
x = np.linspace(xmin, xmax, nx)
y = np.linspace(ymin, ymax, ny)
X, Y = np.meshgrid(x, y)
return X, Y
# 流れのシミュレーション（簡単な例）
def simulate_flow(X, Y, u_max):
U = np.zeros_like(X)
V = np.zeros_like(Y)
# 流れの速度（例として一様な速度場を仮定）
U[:, :] = u_max
return U, V
# メッシュの設定
xmin, xmax = 0, 10
ymin, ymax = 0, 10
nx, ny = 100, 100
X, Y = generate_mesh(xmin, xmax, ymin, ymax, nx, ny)
# 流れのシミュレーション
u_max = 1.0 # 最大流速
U, V = simulate_flow(X, Y, u_max)
# 流れの可視化
plt.figure(figsize=(12, 6))
# 流線のプロット
plt.subplot(1, 2, 1)
plt.title('Streamlines')
plt.streamplot(X, Y, U, V, color='b', linewidth=1.5)
plt.xlabel('X')
plt.ylabel('Y')
# ベクトル場のプロット
plt.subplot(1, 2, 2)
plt.title('Velocity Field')
plt.quiver(X, Y, U, V, scale=20)
plt.xlabel('X')
plt.ylabel('Y')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 風洞実験のシミュレーション（簡単な例）
def airflow_simulation(x, y, velocity):
X, Y = np.meshgrid(x, y)
U = velocity * np.ones_like(X)
V = np.zeros_like(Y)
return U, V
# 水槽実験のシミュレーション（簡単な例）
def water_tank_simulation(x, y, velocity, vortex_strength):
X, Y = np.meshgrid(x, y)
U = velocity * np.ones_like(X) - vortex_strength * (Y - y.mean())
V = vortex_strength * (X - x.mean())
return U, V
# 流速計のデータシミュレーション
def flow_speed_measurement(x, velocity):
return velocity * np.ones_like(x)
# 圧力センサーのデータシミュレーション
def pressure_measurement(x, pressure_gradient):
return pressure_gradient * x
# メッシュの設定
x = np.linspace(0, 10, 100)
y = np.linspace(0, 10, 100)
# 風洞実験のシミュレーション
velocity = 1.0 # 流速 (m/s)
U_airflow, V_airflow = airflow_simulation(x, y, velocity)
# 水槽実験のシミュレーション
vortex_strength = 0.5
U_water_tank, V_water_tank = water_tank_simulation(x, y, velocity, vortex_strength)
# 流速計のデータ
x_flow_speed = np.linspace(0, 10, 100)
flow_speed = flow_speed_measurement(x_flow_speed, velocity)
# 圧力センサーのデータ
pressure_gradient = 0.1
x_pressure = np.linspace(0, 10, 100)
pressure = pressure_measurement(x_pressure, pressure_gradient)
# プロット
plt.figure(figsize=(14, 12))
# 風洞実験のプロット
plt.subplot(2, 2, 1)
plt.title('Airflow Simulation')
plt.streamplot(x, y, U_airflow, V_airflow, color='b')
plt.xlabel('X')
plt.ylabel('Y')
# 水槽実験のプロット
plt.subplot(2, 2, 2)
plt.title('Water Tank Simulation')
plt.streamplot(x, y, U_water_tank, V_water_tank, color='r')
plt.xlabel('X')
plt.ylabel('Y')
# 流速計のデータプロット
plt.subplot(2, 2, 3)
plt.title('Flow Speed Measurement')
plt.plot(x_flow_speed, flow_speed, label='Flow Speed')
plt.xlabel('Position (m)')
plt.ylabel('Speed (m/s)')
plt.legend()
# 圧力センサーのデータプロット
plt.subplot(2, 2, 4)
plt.title('Pressure Measurement')
plt.plot(x_pressure, pressure, label='Pressure')
plt.xlabel('Position (m)')
plt.ylabel('Pressure (Pa)')
plt.legend()
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 降水量と流出量の関係を計算する関数
def runoff_calculation(precipitation, runoff_coefficient):
return precipitation * runoff_coefficient
# データの設定
days = np.arange(1, 31) # 30日間
precipitation = np.random.uniform(0, 10, len(days)) # 0〜10mmのランダムな降水量
runoff_coefficient = 0.6 # 流出係数（例: 0.6）
# 流出量の計算
runoff = runoff_calculation(precipitation, runoff_coefficient)
# プロット
plt.figure(figsize=(12, 6))
# 降水量のプロット
plt.subplot(1, 2, 1)
plt.title('Precipitation')
plt.plot(days, precipitation, label='Precipitation (mm)', color='b')
plt.xlabel('Day')
plt.ylabel('Precipitation (mm)')
plt.legend()
# 流出量のプロット
plt.subplot(1, 2, 2)
plt.title('Runoff')
plt.plot(days, runoff, label='Runoff (mm)', color='r')
plt.xlabel('Day')
plt.ylabel('Runoff (mm)')
plt.legend()
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 擬似塑性流体のパワー低減モデル
def power_law_model(shear_rate, k, n):
"""
shear_rate: 剪断速度
k: 一定
n: パワー則指数
"""
return k * shear_rate**(n-1)
# データの設定
shear_rate = np.linspace(0.01, 10, 100) # 剪断速度の範囲
k = 0.5 # 一定（Pa·s^n）
n = 0.7 # パワー則指数（n < 1は擬似塑性流体）
# 粘度の計算
viscosity = power_law_model(shear_rate, k, n)
# プロット
plt.figure(figsize=(8, 6))
plt.plot(shear_rate, viscosity, label='Power-Law Model', color='b')
plt.xlabel('Shear Rate (1/s)')
plt.ylabel('Viscosity (Pa·s)')
plt.title('Power-Law Model for Non-Newtonian Fluids')
plt.legend()
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 液-気体多相流のシミュレーション
def liquid_gas_flow_simulation(depth, bubble_velocity, liquid_velocity):
"""
depth: 深さ (m)
bubble_velocity: 気泡の上昇速度 (m/s)
liquid_velocity: 液体の速度 (m/s)
"""
time = np.linspace(0, depth / bubble_velocity, 100) # 時間の範囲
bubble_position = bubble_velocity * time
liquid_position = liquid_velocity * time
return time, bubble_position, liquid_position
# データの設定
depth = 10 # 深さ (m)
bubble_velocity = 0.1 # 気泡の上昇速度 (m/s)
liquid_velocity = 0.05 # 液体の速度 (m/s)
# シミュレーション
time, bubble_position, liquid_position = liquid_gas_flow_simulation(depth, bubble_velocity, liquid_velocity)
# プロット
plt.figure(figsize=(8, 6))
plt.plot(time, bubble_position, label='Bubble Position', color='b')
plt.plot(time, liquid_position, label='Liquid Position', color='r')
plt.xlabel('Time (s)')
plt.ylabel('Position (m)')
plt.title('Liquid-Gas Multiphase Flow Simulation')
plt.legend()
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
# データの生成
np.random.seed(0)
velocity = np.random.uniform(0, 5, 100) # 流速 (m/s)
temperature = np.random.uniform(20, 30, 100) # 温度 (°C)
flow_rate = 0.5 * velocity + 0.3 * temperature + np.random.normal(0, 0.5, 100) # 流量 (m³/s)
# モデルのトレーニング
X = np.column_stack((velocity, temperature))
y = flow_rate
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)
model = LinearRegression()
model.fit(X_train, y_train)
# 予測とプロット
y_pred = model.predict(X_test)
plt.scatter(y_test, y_pred, color='blue', label='Predicted vs Actual')
plt.xlabel('Actual Flow Rate')
plt.ylabel('Predicted Flow Rate')
plt.title('Machine Learning Model for Fluid Flow Prediction')
plt.legend()
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Parameters
N = 100 # Grid size
x = np.linspace(-2, 2, N)
y = np.linspace(-2, 2, N)
X, Y = np.meshgrid(x, y)
# Potential function (example: flow around a circular path)
def potential_function(X, Y):
return np.arctan2(Y, X) # Example of complex potential function
# Compute the potential function
phi = potential_function(X, Y)
# Compute the velocity field
U, V = np.gradient(-phi)
# Plot
plt.figure(figsize=(10, 8))
plt.contourf(X, Y, phi, levels=50, cmap='viridis')
plt.colorbar(label='Potential Function')
plt.streamplot(X, Y, U, V, color='white')
plt.title('Potential Function and Velocity Field')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Parameters
N = 100 # Grid size
x = np.linspace(-2, 2, N)
y = np.linspace(-2, 2, N)
X, Y = np.meshgrid(x, y)
# Complex potential function (example: unit circle flow)
def complex_potential(X, Y):
Z = X + 1j * Y
return 1 / Z # Example of complex potential function
# Compute the complex potential function
W = complex_potential(X, Y)
U = np.real(W)
V = np.imag(W)
# Plot
plt.figure(figsize=(10, 8))
plt.contourf(X, Y, np.abs(W), levels=50, cmap='viridis')
plt.colorbar(label='Magnitude of Complex Potential Function')
plt.streamplot(X, Y, U, V, color='white')
plt.title('Complex Potential Function and Velocity Field')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定常流れの速度場
def steady_flow(x, y):
U = np.ones_like(x) # x方向の速度
V = np.zeros_like(y) # y方向の速度
return U, V
# 非定常流れの速度場（時間依存）
def unsteady_flow(x, y, t):
U = np.sin(x) * np.cos(t) # 時間依存の速度
V = np.cos(y) * np.sin(t) # 時間依存の速度
return U, V
# グリッド作成
x = np.linspace(0, 2*np.pi, 100)
y = np.linspace(0, 2*np.pi, 100)
X, Y = np.meshgrid(x, y)
# 定常流れのプロット
U, V = steady_flow(X, Y)
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.quiver(X, Y, U, V)
plt.title('Steady Flow')
plt.xlabel('x')
plt.ylabel('y')
# 非定常流れのプロット
t = 1 # 時間t=1の例
U, V = unsteady_flow(X, Y, t)
plt.subplot(1, 2, 2)
plt.quiver(X, Y, U, V)
plt.title('Unsteady Flow')
plt.xlabel('x')
plt.ylabel('y')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 層流の速度場
def laminar_flow(x, y):
U = np.ones_like(x) # x方向の速度
V = np.zeros_like(y) # y方向の速度
return U, V
# 乱流の速度場
def turbulent_flow(x, y):
np.random.seed(0)
U = np.sin(x) + 0.1 * np.random.randn(*x.shape) # ランダムな擾乱
V = np.cos(y) + 0.1 * np.random.randn(*y.shape) # ランダムな擾乱
return U, V
# グリッド作成
x = np.linspace(0, 2*np.pi, 100)
y = np.linspace(0, 2*np.pi, 100)
X, Y = np.meshgrid(x, y)
# 層流のプロット
U, V = laminar_flow(X, Y)
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.quiver(X, Y, U, V)
plt.title('Laminar Flow')
plt.xlabel('x')
plt.ylabel('y')
# 乱流のプロット
U, V = turbulent_flow(X, Y)
plt.subplot(1, 2, 2)
plt.quiver(X, Y, U, V)
plt.title('Turbulent Flow')
plt.xlabel('x')
plt.ylabel('y')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# 流速場
def velocity_field(XY, t=0):
x, y = XY
U = np.ones_like(x) # x方向の速度
V = np.zeros_like(y) # y方向の速度
return np.array([U, V])
# 流線のプロット
def streamlines(x, y):
def dx_dt(XY, t=0):
x, y = XY
U, V = velocity_field([x, y], t)
return [U, V]
start_points = np.array([[0.5, 0.5], [1.0, 1.0], [1.5, 1.5]])
plt.figure(figsize=(10, 5))
plt.streamplot(X, Y, *velocity_field([X, Y]), color='b')
plt.scatter(start_points[:,0], start_points[:,1], c='r', marker='o')
plt.title('Streamlines')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
# 圧力場のプロット
def pressure_field(x, y):
return np.exp(-x**2 - y**2) # 単純な圧力場
# グリッド作成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 流線と等圧線のプロット
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
streamlines(X, Y)
plt.subplot(1, 2, 2)
plt.contour(X, Y, pressure_field(X, Y))
plt.title('Pressure Contours')
plt.xlabel('x')
plt.ylabel('y')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定義域
x = np.linspace(-2, 2, 400)
y = np.linspace(-2, 2, 400)
X, Y = np.meshgrid(x, y)
# 複素ポテンシャルの定義
def complex_potential_source(z, strength):
return strength * np.log(z)
def complex_potential_vortex(z, strength):
return -strength * np.log(z)
# 複素数平面
Z = X + 1j * Y
# 線源と渦の強度
source_strength = 1.0
vortex_strength = 1.0
# 複素ポテンシャル
W_source = complex_potential_source(Z, source_strength)
W_vortex = complex_potential_vortex(Z, vortex_strength)
# 結果のプロット
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.contourf(X, Y, np.real(W_source), 20, cmap='coolwarm')
plt.colorbar()
plt.title('Real part of Source Potential')
plt.subplot(1, 2, 2)
plt.contourf(X, Y, np.real(W_vortex), 20, cmap='coolwarm')
plt.colorbar()
plt.title('Real part of Vortex Potential')
plt.show()
# 線源、渦、鏡像源の定義
def source(z, strength, z0):
return strength * np.log(z - z0)
def vortex(z, strength, z0):
return -strength * np.log(z - z0)
def mirror_source(z, strength, z0, z_mirror):
return source(z, strength, z0) + source(z, -strength, z_mirror)
def mirror_vortex(z, strength, z0, z_mirror):
return vortex(z, strength, z0) + vortex(z, -strength, z_mirror)
# 鏡像法の設定
source_position = 0.5 + 0.5j
mirror_position = -0.5 - 0.5j
vortex_position = -0.5 + 0.5j
# 線源と渦の強度
source_strength = 1.0
vortex_strength = 1.0
# 複素ポテンシャルの計算
W_combined_source = mirror_source(Z, source_strength, source_position, mirror_position)
W_combined_vortex = mirror_vortex(Z, vortex_strength, vortex_position, mirror_position)
# 結果のプロット
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.contourf(X, Y, np.real(W_combined_source), 20, cmap='coolwarm')
plt.colorbar()
plt.title('Real part of Combined Source Potential')
plt.subplot(1, 2, 2)
plt.contourf(X, Y, np.real(W_combined_vortex), 20, cmap='coolwarm')
plt.colorbar()
plt.title('Real part of Combined Vortex Potential')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定義域
x = np.linspace(-2, 2, 400)
y = np.linspace(-2, 2, 400)
X, Y = np.meshgrid(x, y)
# 円周流れとカースティーン流れの計算
def potential_flow_cylinder(r, gamma):
return -np.log(r)
def potential_flow_kaur(x, y, R, gamma):
r = np.sqrt(x**2 + y**2)
return -np.log(np.sqrt((r**2 + R**2 - 2*R*r*np.cos(gamma)) / (2*R*r)))
# 円の半径と渦の強度
R = 1.0
gamma = 1.0
# 円周流れのポテンシャル計算
r = np.sqrt(X**2 + Y**2)
phi_cylinder = potential_flow_cylinder(r, gamma)
phi_kaur = potential_flow_kaur(X, Y, R, gamma)
# 結果のプロット
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.contourf(X, Y, phi_cylinder, 20, cmap='coolwarm')
plt.colorbar()
plt.title('Potential Flow around Cylinder')
plt.subplot(1, 2, 2)
plt.contourf(X, Y, phi_kaur, 20, cmap='coolwarm')
plt.colorbar()
plt.title('Kaur Flow')
plt.show()
def complex_potential_cylinder(z, R, gamma):
return -np.log(z) + gamma / (2 * np.pi) * np.arctan2(np.imag(z), np.real(z))
# 円柱の半径と渦の強度
R = 1.0
gamma = 1.0
# 複素平面
Z = X + 1j * Y
# 円柱周りのポテンシャル計算
W_cylinder = complex_potential_cylinder(Z, R, gamma)
# 結果のプロット
plt.figure(figsize=(8, 6))
plt.contourf(X, Y, np.real(W_cylinder), 20, cmap='coolwarm')
plt.colorbar()
plt.title('Complex Potential around Cylinder')
plt.show()
def potential_flow_flat_plate(x, y, U_inf):
return U_inf * y
# 無限平板の流れの強度
U_inf = 1.0
# 平板流れのポテンシャル計算
phi_flat_plate = potential_flow_flat_plate(X, Y, U_inf)
# 結果のプロット
plt.figure(figsize=(8, 6))
plt.contourf(X, Y, phi_flat_plate, 20, cmap='coolwarm')
plt.colorbar()
plt.title('Potential Flow over Flat Plate')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定数とパラメータ
nx, ny = 100, 100 # グリッドサイズ
niters = 1000 # イテレーション回数
rho = 1.0 # 密度
nu = 0.1 # 動粘性係数
# 初期条件
u = np.zeros((nx, ny)) # x方向の速度
v = np.zeros((nx, ny)) # y方向の速度
p = np.ones((nx, ny)) * 100 # 圧力
b = np.zeros((nx, ny)) # 中間変数
# 外力
fx = np.zeros((nx, ny))
fy = np.zeros((nx, ny))
# ナビエ-ストークス方程式の数値解法
def solve_ns(u, v, p, b, fx, fy):
un = np.copy(u)
vn = np.copy(v)
pn = np.copy(p)
# 圧力ポアソン方程式
b[1:-1, 1:-1] = (rho * (1 / dt * ((u[2:, 1:-1] - u[:-2, 1:-1]) / 2.0 + (v[1:-1, 2:] - v[1:-1, :-2]) / 2.0) -
((u[2:, 1:-1] - u[:-2, 1:-1]) / 2.0) ** 2 -
2 * ((u[1:-1, 2:] - u[1:-1, :-2]) / 2.0) * ((v[2:, 1:-1] - v[:-2, 1:-1]) / 2.0) -
((v[1:-1, 2:] - v[1:-1, :-2]) / 2.0) ** 2)
p = np.linalg.solve(A, b)
# 速度の更新
u[1:-1, 1:-1] = (un[1:-1, 1:-1] - un[1:-1, 1:-1] * dt / dx * (un[1:-1, 1:-1] - un[:-2, 1:-1]) -
vn[1:-1, 1:-1] * dt / dy * (un[1:-1, 1:-1] - un[1:-1, :-2]) -
dt / (2 * dx * rho) * (p[2:, 1:-1] - p[:-2, 1:-1]) +
nu * ((un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[:-2, 1:-1]) / dx**2 +
(un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, :-2]) / dy**2))
v[1:-1, 1:-1] = (vn[1:-1, 1:-1] - un[1:-1, 1:-1] * dt / dx * (vn[1:-1, 1:-1] - vn[:-2, 1:-1]) -
vn[1:-1, 1:-1] * dt / dy * (vn[1:-1, 1:-1] - vn[1:-1, :-2]) -
dt / (2 * dy * rho) * (p[1:-1, 2:] - p[1:-1, :-2]) +
nu * ((vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[:-2, 1:-1]) / dx**2 +
(vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, :-2]) / dy**2))
return u, v, p
# イテレーションのループ
for _ in range(niters):
u, v, p = solve_ns(u, v, p, b, fx, fy)
# 結果の表示
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.title('Velocity U')
plt.imshow(u, origin='lower', cmap='jet')
plt.colorbar()
plt.subplot(1, 2, 2)
plt.title('Pressure')
plt.imshow(p, origin='lower', cmap='jet')
plt.colorbar()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve
# 定数とパラメータ
nx, ny = 100, 100 # グリッドサイズ
niters = 100 # イテレーション回数
dx, dy = 1.0, 1.0 # 空間間隔
dt = 0.1 # 時間間隔
rho = 1.0 # 密度
nu = 0.1 # 動粘性係数
# 初期条件
u = np.zeros((nx, ny)) # x方向の速度
v = np.zeros((nx, ny)) # y方向の速度
p = np.ones((nx, ny)) * 100 # 圧力
b = np.zeros((nx, ny)) # 中間変数
# 外力
fx = np.zeros((nx, ny))
fy = np.zeros((nx, ny))
# ナビエ-ストークス方程式の数値解法
def solve_ns(u, v, p, b, fx, fy):
un = np.copy(u)
vn = np.copy(v)
pn = np.copy(p)
# 圧力ポアソン方程式（簡略化）
b[1:-1, 1:-1] = (rho * (1 / dt * ((u[2:, 1:-1] - u[:-2, 1:-1]) / 2.0 + (v[1:-1, 2:] - v[1:-1, :-2]) / 2.0) -
((u[2:, 1:-1] - u[:-2, 1:-1]) / 2.0) ** 2 -
2 * ((u[1:-1, 2:] - u[1:-1, :-2]) / 2.0) * ((v[2:, 1:-1] - v[:-2, 1:-1]) / 2.0) -
((v[1:-1, 2:] - v[1:-1, :-2]) / 2.0) ** 2)
A = np.eye(nx * ny) # 単位行列（仮のもの）
p = np.linalg.solve(A, b.flatten()).reshape((nx, ny))
# 速度の更新
u[1:-1, 1:-1] = (un[1:-1, 1:-1] - un[1:-1, 1:-1] * dt / dx * (un[1:-1, 1:-1] - un[:-2, 1:-1]) -
vn[1:-1, 1:-1] * dt / dy * (un[1:-1, 1:-1] - un[1:-1, :-2]) -
dt / (2 * dx * rho) * (p[2:, 1:-1] - p[:-2, 1:-1]) +
nu * ((un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[:-2, 1:-1]) / dx**2 +
(un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, :-2]) / dy**2))
v[1:-1, 1:-1] = (vn[1:-1, 1:-1] - un[1:-1, 1:-1] * dt / dx * (vn[1:-1, 1:-1] - vn[:-2, 1:-1]) -
vn[1:-1, 1:-1] * dt / dy * (vn[1:-1, 1:-1] - vn[1:-1, :-2]) -
dt / (2 * dy * rho) * (p[1:-1, 2:] - p[1:-1, :-2]) +
nu * ((vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[:-2, 1:-1]) / dx**2 +
(vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, :-2]) / dy**2))
return u, v, p
# イテレーションのループ
for _ in range(niters):
u, v, p = solve_ns(u, v, p, b, fx, fy)
# プロット
plt.figure(figsize=(14, 12))
# 1. 初期速度場
plt.subplot(2, 2, 1)
plt.title('Initial Velocity Field')
plt.quiver(u, v, scale=5, scale_units='inches', angles='xy')
plt.xlabel('x')
plt.ylabel('y')
# 2. 速度場の分布（速度の大きさ）
plt.subplot(2, 2, 2)
plt.title('Velocity Magnitude')
speed = np.sqrt(u**2 + v**2)
plt.imshow(speed, origin='lower', cmap='viridis')
plt.colorbar(label='Speed')
plt.xlabel('x')
plt.ylabel('y')
# 3. 圧力場の分布
plt.subplot(2, 2, 3)
plt.title('Pressure Field')
plt.imshow(p, origin='lower', cmap='plasma')
plt.colorbar(label='Pressure')
plt.xlabel('x')
plt.ylabel('y')
# 4. 速度場の時間経過による変化（例として1回目と最終回の比較）
plt.subplot(2, 2, 4)
plt.title('Velocity Field Evolution')
plt.quiver(u, v, scale=5, scale_units='inches', angles='xy')
plt.xlabel('x')
plt.ylabel('y')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定数とパラメータ
nx, ny = 100, 100 # グリッドサイズ
niters = 100 # イテレーション回数
dx, dy = 1.0, 1.0 # 空間間隔
dt = 0.1 # 時間間隔
rho = 1.0 # 密度
nu = 0.1 # 動粘性係数
# 初期条件
u = np.zeros((nx, ny), dtype=complex) # x方向の速度（複素数）
v = np.zeros((nx, ny), dtype=complex) # y方向の速度（複素数）
p = np.ones((nx, ny)) * 100 # 圧力
b = np.zeros((nx, ny)) # 中間変数
# 外力
fx = np.zeros((nx, ny))
fy = np.zeros((nx, ny))
# ナビエ-ストークス方程式の数値解法
def solve_ns(u, v, p, b, fx, fy):
un = np.copy(u)
vn = np.copy(v)
pn = np.copy(p)
# 圧力ポアソン方程式（簡略化）
b[1:-1, 1:-1] = (rho * (1 / dt * ((u[2:, 1:-1] - u[:-2, 1:-1]) / 2.0 + (v[1:-1, 2:] - v[1:-1, :-2]) / 2.0) -
((u[2:, 1:-1] - u[:-2, 1:-1]) / 2.0) ** 2 -
2 * ((u[1:-1, 2:] - u[1:-1, :-2]) / 2.0) * ((v[2:, 1:-1] - v[:-2, 1:-1]) / 2.0) -
((v[1:-1, 2:] - v[1:-1, :-2]) / 2.0) ** 2)
A = np.eye(nx * ny) # 単位行列（仮のもの）
p = np.linalg.solve(A, b.flatten()).reshape((nx, ny))
# 速度の更新
u[1:-1, 1:-1] = (un[1:-1, 1:-1] - un[1:-1, 1:-1] * dt / dx * (un[1:-1, 1:-1] - un[:-2, 1:-1]) -
vn[1:-1, 1:-1] * dt / dy * (un[1:-1, 1:-1] - un[1:-1, :-2]) -
dt / (2 * dx * rho) * (p[2:, 1:-1] - p[:-2, 1:-1]) +
nu * ((un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[:-2, 1:-1]) / dx**2 +
(un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, :-2]) / dy**2))
v[1:-1, 1:-1] = (vn[1:-1, 1:-1] - un[1:-1, 1:-1] * dt / dx * (vn[1:-1, 1:-1] - vn[:-2, 1:-1]) -
vn[1:-1, 1:-1] * dt / dy * (vn[1:-1, 1:-1] - vn[1:-1, :-2]) -
dt / (2 * dy * rho) * (p[1:-1, 2:] - p[1:-1, :-2]) +
nu * ((vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[:-2, 1:-1]) / dx**2 +
(vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, :-2]) / dy**2))
return u, v, p
# イテレーションのループ
for _ in range(niters):
u, v, p = solve_ns(u, v, p, b, fx, fy)
# プロット
plt.figure(figsize=(14, 12))
# 1. 複素数平面での速度場の可視化
plt.subplot(2, 2, 1)
plt.title('Complex Velocity Field')
plt.quiver(np.real(u), np.imag(u), scale=5, scale_units='inches', angles='xy')
plt.xlabel('x')
plt.ylabel('y')
# 2. 速度場の大きさのヒートマップ
plt.subplot(2, 2, 2)
plt.title('Velocity Magnitude')
speed = np.abs(u + 1j * v) # 複素数としての速度場の大きさ
plt.imshow(speed, origin='lower', cmap='viridis')
plt.colorbar(label='Speed')
plt.xlabel('x')
plt.ylabel('y')
# 3. 複素数平面上の圧力場の可視化
plt.subplot(2, 2, 3)
plt.title('Pressure Field')
plt.imshow(p, origin='lower', cmap='plasma')
plt.colorbar(label='Pressure')
plt.xlabel('x')
plt.ylabel('y')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Parameters
nx, ny = 41, 41
nt = 500
nit = 50
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
rho = 1
nu = .1
dt = .001
# Initialize variables
u = np.zeros((ny, nx))
v = np.zeros((ny, nx))
p = np.zeros((ny, nx))
b = np.zeros((ny, nx))
# Helper functions
def build_up_b(b, rho, dt, u, v, dx, dy):
b[1:-1, 1:-1] = (rho * (1 / dt *
((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx) +
(v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) -
((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx))**2 -
2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dy) *
(v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dx))-
((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy))**2))
return b
def pressure_poisson(p, dx, dy, b):
pn = np.empty_like(p)
pn = p.copy()
for q in range(nit):
pn = p.copy()
p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +
(pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2) /
(2 * (dx**2 + dy**2)) -
dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * b[1:-1, 1:-1])
p[:, -1] = p[:, -2] # dp/dx = 0 at x = 2
p[0, :] = p[1, :] # dp/dy = 0 at y = 0
p[:, 0] = p[:, 1] # dp/dx = 0 at x = 0
p[-1, :] = 0 # p = 0 at y = 2
return p
def cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu):
un = np.empty_like(u)
vn = np.empty_like(v)
b = np.zeros((ny, nx))
for n in range(nt):
un = u.copy()
vn = v.copy()
b = build_up_b(b, rho, dt, u, v, dx, dy)
p = pressure_poisson(p, dx, dy, b)
u[1:-1, 1:-1] = (un[1:-1, 1:-1] -
un[1:-1, 1:-1] * dt / dx *
(un[1:-1, 1:-1] - un[1:-1, 0:-2]) -
vn[1:-1, 1:-1] * dt / dy *
(un[1:-1, 1:-1] - un[0:-2, 1:-1]) -
dt / (2 * rho * dx) * (p[1:-1, 2:] - p[1:-1, 0:-2]) +
nu * (dt / dx**2 *
(un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
dt / dy**2 *
(un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])))
v[1:-1, 1:-1] = (vn[1:-1, 1:-1] -
un[1:-1, 1:-1] * dt / dx *
(vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -
vn[1:-1, 1:-1] * dt / dy *
(vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) -
dt / (2 * rho * dy) * (p[2:, 1:-1] - p[0:-2, 1:-1]) +
nu * (dt / dx**2 *
(vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
dt / dy**2 *
(vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))
u[0, :] = 0
u[:, 0] = 0
u[:, -1] = 0
u[-1, :] = 1 # set velocity on cavity lid equal to 1
v[0, :] = 0
v[:, 0] = 0
v[:, -1] = 0
v[-1, :] = 0
return u, v, p
# Running the simulation
u, v, p = cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu)
# Plotting
plt.figure(figsize=(11, 7), dpi=100)
plt.contourf(np.linspace(0, 2, nx), np.linspace(0, 2, ny), p, alpha=0.5, cmap=plt.cm.viridis)
plt.colorbar()
plt.contour(np.linspace(0, 2, nx), np.linspace(0, 2, ny), p, cmap=plt.cm.viridis)
plt.quiver(np.linspace(0, 2, nx), np.linspace(0, 2, ny), u, v)
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Cavity Flow with Navier-Stokes Equations')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Grid parameters
nx, ny = 50, 50
dx, dy = 1.0 / nx, 1.0 / ny
dt = 0.001
nu = 0.01
rho = 1.0
# Initialize fields
u = np.zeros((nx, ny))
v = np.zeros((nx, ny))
k = np.ones((nx, ny)) * 0.01
epsilon = np.ones((nx, ny)) * 0.01
nut = np.ones((nx, ny)) * 0.01
# Simulation parameters
Cmu = 0.09
sigma_k = 1.0
sigma_epsilon = 1.3
C1 = 1.44
C2 = 1.92
# Time-stepping loop
for n in range(1000):
# Update turbulent viscosity
nut = Cmu * k**2 / epsilon
# Update k and epsilon using finite difference
k[1:-1, 1:-1] += dt * (
nut[1:-1, 1:-1] * (np.gradient(np.gradient(k, dx, axis=0), dx, axis=0) + np.gradient(np.gradient(k, dy, axis=1), dy, axis=1)) -
u[1:-1, 1:-1] * np.gradient(k, dx, axis=0) -
v[1:-1, 1:-1] * np.gradient(k, dy, axis=1) +
nut[1:-1, 1:-1] * (np.gradient(np.gradient(k, dx, axis=0), dx, axis=0) + np.gradient(np.gradient(k, dy, axis=1), dy, axis=1)) -
epsilon[1:-1, 1:-1]
)
epsilon[1:-1, 1:-1] += dt * (
nut[1:-1, 1:-1] * (np.gradient(np.gradient(epsilon, dx, axis=0), dx, axis=0) + np.gradient(np.gradient(epsilon, dy, axis=1), dy, axis=1)) -
u[1:-1, 1:-1] * np.gradient(epsilon, dx, axis=0) -
v[1:-1, 1:-1] * np.gradient(epsilon, dy, axis=1) +
C1 * epsilon[1:-1, 1:-1] / k[1:-1, 1:-1] * nut[1:-1, 1:-1] * (np.gradient(u[1:-1, 1:-1], dx, axis=0) + np.gradient(v[1:-1, 1:-1], dy, axis=1)) -
C2 * epsilon[1:-1, 1:-1]**2 / k[1:-1, 1:-1]
)
# Plotting
plt.figure()
plt.contourf(k.T, cmap='jet')
plt.colorbar()
plt.title('Turbulent Kinetic Energy (k)')
plt.show()
plt.figure()
plt.contourf(epsilon.T, cmap='jet')
plt.colorbar()
plt.title('Dissipation Rate (epsilon)')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# Grid parameters
nx, ny = 50, 50
dx, dy = 1.0 / nx, 1.0 / ny
dt = 0.001
rho_l = 1000.0 # Density of liquid
rho_g = 1.0 # Density of gas
# Initialize fields
u = np.zeros((nx, ny))
v = np.zeros((nx, ny))
alpha_l = np.zeros((nx, ny))
alpha_g = np.zeros((nx, ny))
alpha_l[:, :] = 1.0 # Initially filled with liquid
# Introduce gas bubble
alpha_g[20:30, 20:30] = 1.0
alpha_l[20:30, 20:30] = 0.0
# Time-stepping loop
for n in range(1000):
# Compute new alpha fields
alpha_l[1:-1, 1:-1] += dt * (
- u[1:-1, 1:-1] * np.gradient(alpha_l, dx, axis=0)[1:-1, 1:-1] -
v[1:-1, 1:-1] * np.gradient(alpha_l, dy, axis=1)[1:-1, 1:-1]
)
alpha_g = 1.0 - alpha_l # Complementary relation
# Update velocity fields (simplified)
u[1:-1, 1:-1] += dt * (-alpha_g[1:-1, 1:-1] * np.gradient(p, dx, axis=0)[1:-1, 1:-1])
v[1:-1, 1:-1] += dt * (-alpha_g[1:-1, 1:-1] * np.gradient(p, dy, axis=1)[1:-1, 1:-1])
# Plotting
plt.figure()
plt.contourf(alpha_l.T, cmap='Blues')
plt.colorbar()
plt.title('Liquid Volume Fraction')
plt.show()
plt.figure()
plt.contourf(alpha_g.T, cmap='Reds')
plt.colorbar()
plt.title('Gas Volume Fraction')
plt.show()
import numpy as np
def hagen_poiseuille(delta_p, mu, L, Q, R):
"""
ハーゲン・ポアズイユの式を使用して圧力損失を計算する関数
delta_p : 圧力損失 (Pa)
mu : 流体の動粘度 (Pa·s)
L : 管の長さ (m)
Q : 体積流量 (m³/s)
R : 管の半径 (m)
"""
return 8 * mu * L * Q / (np.pi * R**4)
# 例としてのパラメータ
mu = 0.001 # 流体の動粘度 (Pa·s)
L = 1.0 # 管の長さ (m)
Q = 0.01 # 体積流量 (m³/s)
R = 0.01 # 管の半径 (m)
# 圧力損失を計算
delta_p = hagen_poiseuille(0, mu, L, Q, R)
print(f"圧力損失: {delta_p:.2f} Pa")
def calculate_lambda(Re):
"""
層流における管摩擦係数を計算する関数
Re : レイノルズ数
"""
return 64 / Re
def reynolds_number(rho, v, D, mu):
"""
レイノルズ数を計算する関数
rho : 流体の密度 (kg/m³)
v : 流速 (m/s)
D : 管の直径 (m)
mu : 流体の動粘度 (Pa·s)
"""
return (rho * v * D) / mu
# 例としてのパラメータ
rho = 1000 # 流体の密度 (kg/m³)
v = 0.5 # 流速 (m/s)
D = 0.02 # 管の直径 (m)
mu = 0.001 # 流体の動粘度 (Pa·s)
# レイノルズ数を計算
Re = reynolds_number(rho, v, D, mu)
# 管摩擦係数を計算
lambda_value = calculate_lambda(Re)
print(f"レイノルズ数: {Re:.2f}")
print(f"管摩擦係数 λ: {lambda_value:.2f}")
import numpy as np
import matplotlib.pyplot as plt
# 流体の性質
density = 1000 # kg/m^3 (水の密度)
viscosity = 0.001 # Pa.s (水の粘度)
# 流体の静力学
def hydrostatic_pressure(depth, density, g=9.81):
return density * g * depth
depths = np.linspace(0, 10, 100)
pressures = hydrostatic_pressure(depths, density)
plt.figure()
plt.plot(depths, pressures)
plt.xlabel('Depth (m)')
plt.ylabel('Pressure (Pa)')
plt.title('Hydrostatic Pressure')
plt.grid(True)
plt.show()
# 流れの力学
def bernoulli(p1, v1, p2, h1, h2, density, g=9.81):
v2 = np.sqrt(2 * (p1 - p2) / density + v1**2 + 2 * g * (h1 - h2))
return v2
p1 = 100000 # Pa
v1 = 0 # m/s
p2 = 50000 # Pa
h1 = 10 # m
h2 = 0 # m
v2 = bernoulli(p1, v1, p2, h1, h2, density)
print(f'Exit velocity: {v2:.2f} m/s')
# 管内の流れ
def laminar_flow_velocity(radius, viscosity, pressure_gradient, length):
return (pressure_gradient * radius**2) / (4 * viscosity * length)
radius = 0.01 # m
length = 1 # m
pressure_gradient = 1000 # Pa/m
velocity = laminar_flow_velocity(radius, viscosity, pressure_gradient, length)
print(f'Laminar flow velocity: {velocity:.2f} m/s')
# バルブ特性
def valve_flow_coefficient(flow_rate, density, pressure_drop):
return flow_rate / np.sqrt(pressure_drop / density)
flow_rate = 0.01 # m^3/s
pressure_drop = 5000 # Pa
Cv = valve_flow_coefficient(flow_rate, density, pressure_drop)
print(f'Valve flow coefficient (Cv): {Cv:.2f}')
# プロット
plt.figure()
plt.plot(depths, hydrostatic_pressure(depths, density), label='Hydrostatic Pressure')
plt.axhline(y=p1, color='r', linestyle='--', label='Initial Pressure')
plt.axhline(y=p2, color='b', linestyle='--', label='Final Pressure')
plt.xlabel('Depth (m)')
plt.ylabel('Pressure (Pa)')
plt.title('Fluid Dynamics Summary')
plt.legend()
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
# 1.1 流体力学と流体の性質
# 流体の密度と粘度に関する基本的な計算
def fluid_properties(density, viscosity, velocity, characteristic_length):
reynolds_number = (density * velocity * characteristic_length) / viscosity
return reynolds_number
# パラメータ設定
density = 1000 # 水の密度 (kg/m^3)
viscosity = 0.001 # 水の動粘度 (Pa.s)
velocity = 2.0 # 流速 (m/s)
characteristic_length = 0.05 # 特徴長 (m)
# レイノルズ数の計算
reynolds_number = fluid_properties(density, viscosity, velocity, characteristic_length)
print(f'Reynolds Number: {reynolds_number}')
# 1.2 流体の圧縮性と表面張力
# 圧縮性の計算 (例: 音速の計算)
def sound_speed(density, bulk_modulus):
return np.sqrt(bulk_modulus / density)
# 表面張力の計算 (例: ライマンの法則に基づく)
def surface_tension(density, gravity, radius):
return 2 * density * gravity * radius
# パラメータ設定
bulk_modulus = 2.2e9 # 水の体積弾性率 (Pa)
gravity = 9.81 # 重力加速度 (m/s^2)
radius = 0.01 # 液滴の半径 (m)
# 音速の計算
speed_of_sound = sound_speed(density, bulk_modulus)
print(f'Speed of Sound: {speed_of_sound} m/s')
# 表面張力の計算
surface_tension_value = surface_tension(density, gravity, radius)
print(f'Surface Tension: {surface_tension_value} N/m')
# 1.3 流れのとらえ方
# 流線の可視化 (簡単な2D流れの例)
def plot_streamlines():
x, y = np.meshgrid(np.linspace(-1, 1, 100), np.linspace(-1, 1, 100))
u = -y # 流れのベクトル場 (例: 回転流)
v = x
plt.figure(figsize=(8, 6))
plt.streamplot(x, y, u, v, density=1.5, color='b')
plt.title('Streamlines of a Simple Flow')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
plot_streamlines()
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
# 2.1 力，応力，圧力
# 圧力の計算 (静水圧)
def pressure(depth, density, gravity):
return density * gravity * depth
# 力の計算 (圧力による力)
def force_on_surface(pressure, area):
return pressure * area
# 応力の計算 (単位面積あたりの力)
def stress(force, area):
return force / area
# パラメータ設定
density = 1000 # 水の密度 (kg/m^3)
gravity = 9.81 # 重力加速度 (m/s^2)
depth = 10 # 深さ (m)
area = 2.0 # 表面積 (m^2)
# 圧力の計算
pressure_value = pressure(depth, density, gravity)
print(f'Pressure at depth {depth} m: {pressure_value} Pa')
# 力の計算
force_value = force_on_surface(pressure_value, area)
print(f'Force on surface: {force_value} N')
# 応力の計算
stress_value = stress(force_value, area)
print(f'Stress on surface: {stress_value} Pa')
# 2.2 マノメータ
# マノメータの差圧計算
def manometer_reading(pressure1, pressure2):
return pressure1 - pressure2
# パラメータ設定
pressure1 = 101325 # 大気圧 (Pa)
pressure2 = 100000 # 測定圧力 (Pa)
# マノメータの読み取り値
manometer_value = manometer_reading(pressure1, pressure2)
print(f'Manometer Reading: {manometer_value} Pa')
# 2.3 全圧力と圧力中心
# 全圧力の計算
def total_pressure(density, gravity, height):
return density * gravity * height
# 圧力中心の計算 (面積分での重心)
def pressure_center(depth, area, density, gravity):
return depth / 2 # 簡略化された計算 (深さの半分)
# パラメータ設定
height = 5 # 浸水深 (m)
# 全圧力の計算
total_pressure_value = total_pressure(density, gravity, height)
print(f'Total Pressure: {total_pressure_value} N/m^2')
# 圧力中心の計算
pressure_center_value = pressure_center(height, area, density, gravity)
print(f'Pressure Center Depth: {pressure_center_value} m')
# 2.4 浮力と浮揚体の安定性
# 浮力の計算 (アルキメデスの原理)
def buoyant_force(volume, density, gravity):
return density * gravity * volume
# 浮揚体の安定性 (トルク計算)
def stability_torque(volume, density, gravity, height_center):
return buoyant_force(volume, density, gravity) * height_center
# パラメータ設定
volume = 1.0 # 浸水体積 (m^3)
height_center = height / 2 # 浮揚体の重心
# 浮力の計算
buoyant_force_value = buoyant_force(volume, density, gravity)
print(f'Buoyant Force: {buoyant_force_value} N')
# 安定性トルクの計算
stability_torque_value = stability_torque(volume, density, gravity, height_center)
print(f'Stability Torque: {stability_torque_value} Nm')
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# 3.1 流れの速度と流れる量
# 流速と流量の計算
def flow_rate(velocity, area):
return velocity * area
# パラメータ設定
velocity = 3.0 # 流速 (m/s)
area = 0.1 # 流れる断面積 (m^2)
# 流量の計算
flow_rate_value = flow_rate(velocity, area)
print(f'Flow Rate: {flow_rate_value} m^3/s')
# 3.2 流れの状態
# 流れの状態を示す例として、連続の方程式を使用
def continuity_equation(velocity1, area1, velocity2, area2):
return velocity1 * area1 == velocity2 * area2
# パラメータ設定
velocity1 = 3.0 # 流速 (m/s) - 断面1
area1 = 0.1 # 断面1の面積 (m^2)
velocity2 = 1.5 # 流速 (m/s) - 断面2
area2 = 0.2 # 断面2の面積 (m^2)
# 連続の方程式の確認
continuity_check = continuity_equation(velocity1, area1, velocity2, area2)
print(f'Continuity Equation Satisfied: {continuity_check}')
# 3.3 一次元流れの場合の基礎方程式
# 一次元流れの例として、ベルヌーイの方程式を使用
def bernoulli_equation(p1, p2, velocity1, velocity2, density, height1, height2):
return p1 + 0.5 * density * velocity1**2 + density * 9.81 * height1 == p2 + 0.5 * density * velocity2**2 + density * 9.81 * height2
# パラメータ設定
p1 = 100000 # 圧力1 (Pa)
p2 = 95000 # 圧力2 (Pa)
density = 1000 # 水の密度 (kg/m^3)
height1 = 5 # 高さ1 (m)
height2 = 2 # 高さ2 (m)
# ベルヌーイ方程式の確認
bernoulli_check = bernoulli_equation(p1, p2, velocity1, velocity2, density, height1, height2)
print(f'Bernoulli Equation Satisfied: {bernoulli_check}')
# 3.4 二次元流れの場合の基礎方程式
# 二次元流れの例として、ナビエ-ストークス方程式を用いた数値計算
def navier_stokes(x, y, u, v, mu, density):
# Navier-Stokes equations are complex and often solved numerically
# Here is a simple example for illustrative purposes
return (density * (u**2 + v**2)) - mu * (u**2 + v**2)
# パラメータ設定
mu = 0.001 # 動粘度 (Pa.s)
# 数値計算のための格子と初期条件設定
x = np.linspace(0, 1, 100)
y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(x, y)
U = np.sin(np.pi * X) * np.cos(np.pi * Y)
V = -np.cos(np.pi * X) * np.sin(np.pi * Y)
# 2D流れの可視化
plt.figure(figsize=(8, 6))
plt.quiver(X, Y, U, V, scale=5)
plt.title('2D Flow Field Visualization')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 4.1 流体におけるエネルギー保存則
def bernoulli_equation(p1, p2, v1, v2, density, h1, h2):
"""
ベルヌーイの方程式を使用して、流体のエネルギー保存を確認する関数。
p1, p2 : 圧力 (Pa)
v1, v2 : 流速 (m/s)
density : 密度 (kg/m^3)
h1, h2 : 高さ (m)
"""
return p1 + 0.5 * density * v1**2 + density * 9.81 * h1 == p2 + 0.5 * density * v2**2 + density * 9.81 * h2
# パラメータ設定
p1 = 100000 # 圧力1 (Pa)
p2 = 95000 # 圧力2 (Pa)
v1 = 3.0 # 流速1 (m/s)
v2 = 1.5 # 流速2 (m/s)
density = 1000 # 水の密度 (kg/m^3)
h1 = 5 # 高さ1 (m)
h2 = 2 # 高さ2 (m)
# ベルヌーイ方程式の確認
bernoulli_check = bernoulli_equation(p1, p2, v1, v2, density, h1, h2)
print(f'Bernoulli Equation Satisfied: {bernoulli_check}')
# 4.2 ベルヌーイの定理の応用
def velocity_from_pressure(p1, p2, density, v1, h1, h2):
"""
ベルヌーイの方程式から流速を計算する関数。
p1, p2 : 圧力 (Pa)
density : 密度 (kg/m^3)
v1 : 流速1 (m/s)
h1, h2 : 高さ (m)
"""
v2_squared = (2 * (p1 - p2) + density * (v1**2 + 2 * 9.81 * (h1 - h2))) / density
return np.sqrt(v2_squared) if v2_squared > 0 else 0
# 流速の計算
v2_calculated = velocity_from_pressure(p1, p2, density, v1, h1, h2)
print(f'Calculated Velocity v2: {v2_calculated:.2f} m/s')
# 4.3 流体の速度・流量の測定
def plot_velocity_pressure_relationship(p1, p2, density, v1, h1, h2):
"""
流体の速度と圧力の関係をプロットする関数。
"""
pressures = np.linspace(p2, p1, 100)
velocities = [velocity_from_pressure(p1, p, density, v1, h1, h2) for p in pressures]
plt.figure(figsize=(10, 6))
plt.plot(pressures, velocities, label='Velocity vs Pressure')
plt.title('Relationship between Pressure and Velocity')
plt.xlabel('Pressure (Pa)')
plt.ylabel('Velocity (m/s)')
plt.legend()
plt.grid(True)
plt.show()
# 圧力と流速の関係をプロット
plot_velocity_pressure_relationship(p1, p2, density, v1, h1, h2)
import numpy as np
# ダルシー・ワイスバッハの式による圧力損失
def darcy_weisbach_loss(flow_rate, diameter, length, roughness, viscosity, density):
area = np.pi * (diameter / 2) ** 2
velocity = flow_rate / area
reynolds_number = (density * velocity * diameter) / viscosity
if reynolds_number < 2000: # 層流の場合
friction_factor = 64 / reynolds_number
else: # 乱流の場合
friction_factor = 0.079 / reynolds_number**0.25
pressure_loss = friction_factor * (length / diameter) * (density * velocity**2 / 2)
return pressure_loss
# 層流の場合の摩擦損失
def laminar_friction_loss(flow_rate, diameter, length, viscosity, density):
return darcy_weisbach_loss(flow_rate, diameter, length, 0, viscosity, density)
# 乱流の場合の摩擦損失
def turbulent_friction_loss(flow_rate, diameter, length, roughness, viscosity, density):
return darcy_weisbach_loss(flow_rate, diameter, length, roughness, viscosity, density)
# 局所損失を考慮した総損失
def local_losses(flow_rate, diameter, length, roughness, viscosity, density, local_loss_coefficient):
friction_loss = darcy_weisbach_loss(flow_rate, diameter, length, roughness, viscosity, density)
return friction_loss + local_loss_coefficient * (density * flow_rate**2 / (2 * diameter))
# 管路全体の総損失
def total_pipeline_loss(flow_rate, diameter, length, roughness, viscosity, density, local_loss_coefficient):
return darcy_weisbach_loss(flow_rate, diameter, length, roughness, viscosity, density) + \
local_loss_coefficient * (density * flow_rate**2 / (2 * diameter))
# 抗力を計算する関数
def drag_force(density, velocity, area, drag_coefficient):
return 0.5 * density * velocity**2 * area * drag_coefficient
# 揚力を計算する関数
def lift_force(density, velocity, area, lift_coefficient):
return 0.5 * density * velocity**2 * area * lift_coefficient
# 粘性流体による流動の圧力損失を計算する関数
def viscous_flow_loss(flow_rate, viscosity, length, diameter):
area = np.pi * (diameter / 2)**2
velocity = flow_rate / area
reynolds_number = (density * velocity * diameter) / viscosity
return (32 * viscosity * length * velocity) / (density * diameter**2)
# パラメータ設定
flow_rate = 0.01 # 流量 (m^3/s)
diameter = 0.05 # 管の直径 (m)
length = 10 # 管の長さ (m)
roughness = 0.0001 # 管の粗さ (m)
viscosity = 1e-6 # 動粘度 (m^2/s)
density = 1000 # 流体の密度 (kg/m^3)
local_loss_coefficient = 0.5 # 局所損失係数
area = 1.0 # 物体の投影面積 (m^2)
drag_coefficient = 0.47 # 抗力係数
lift_coefficient = 1.2 # 揚力係数
# 圧力損失の計算
pressure_loss = darcy_weisbach_loss(flow_rate, diameter, length, roughness, viscosity, density)
print(f'Pressure Loss: {pressure_loss:.2f} Pa')
# 層流の摩擦損失の計算
laminar_loss = laminar_friction_loss(flow_rate, diameter, length, viscosity, density)
print(f'Laminar Friction Loss: {laminar_loss:.2f} Pa')
# 乱流の摩擦損失の計算
turbulent_loss = turbulent_friction_loss(flow_rate, diameter, length, roughness, viscosity, density)
print(f'Turbulent Friction Loss: {turbulent_loss:.2f} Pa')
# 局所損失を考慮した総損失の計算
total_loss = local_losses(flow_rate, diameter, length, roughness, viscosity, density, local_loss_coefficient)
print(f'Total Loss with Local Losses: {total_loss:.2f} Pa')
# 総損失の計算
pipeline_loss = total_pipeline_loss(flow_rate, diameter, length, roughness, viscosity, density, local_loss_coefficient)
print(f'Total Pipeline Loss: {pipeline_loss:.2f} Pa')
# 抗力と揚力の計算
drag = drag_force(density, 30, area, drag_coefficient)
lift = lift_force(density, 30, area, lift_coefficient)
print(f'Drag Force: {drag:.2f} N')
print(f'Lift Force: {lift:.2f} N')
# 粘性流体による流動の圧力損失の計算
viscous_loss = viscous_flow_loss(flow_rate, viscosity, length, diameter)
print(f'Viscous Flow Loss: {viscous_loss:.2f} Pa')
import numpy as np
import matplotlib.pyplot as plt
# ポテンシャル関数の定義（例: 2Dポテンシャル流れ）
def potential_function(x, y):
return x**2 - y**2
# 流れの場の計算
def velocity_field(x, y):
dx, dy = np.gradient(potential_function(x, y))
return -dx, -dy
# グリッドの作成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# ポテンシャルと速度場の計算
phi = potential_function(X, Y)
U, V = velocity_field(X, Y)
# ポテンシャルのプロット
plt.figure(figsize=(12, 6))
# ポテンシャルの等高線
plt.contour(X, Y, phi, levels=20, cmap='jet')
plt.colorbar(label='Potential Function')
# 流れの場のベクトルプロット
plt.quiver(X, Y, U, V, color='black', scale=20)
plt.title('Potential Flow and Velocity Field')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 速度ポテンシャルと流れ関数の定義
def velocity_potential(x, y):
return x**2 - y**2
def stream_function(x, y):
return x * y
# グリッドの作成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# ポテンシャルと流れ関数の計算
Phi = velocity_potential(X, Y)
Psi = stream_function(X, Y)
# プロット
plt.figure(figsize=(12, 6))
# 速度ポテンシャルの等高線
plt.contour(X, Y, Phi, levels=20, cmap='jet')
plt.colorbar(label='Velocity Potential (Phi)')
# 流れ関数の等高線
plt.contour(X, Y, Psi, levels=20, colors='black')
plt.colorbar(label='Stream Function (Psi)')
plt.title('Velocity Potential and Stream Function')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 複素速度ポテンシャルの定義
def complex_velocity_potential(x, y):
phi = x**2 - y**2
psi = x * y
return phi + 1j * psi
# グリッドの作成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 複素速度ポテンシャルの計算
Z = complex_velocity_potential(X, Y)
Phi_real = np.real(Z)
Psi_imag = np.imag(Z)
# プロット
plt.figure(figsize=(12, 6))
# 実部の等高線
plt.contour(X, Y, Phi_real, levels=20, cmap='jet')
plt.colorbar(label='Real Part of Complex Potential')
# 虚部の等高線
plt.contour(X, Y, Psi_imag, levels=20, colors='black')
plt.colorbar(label='Imaginary Part of Complex Potential')
plt.title('Complex Velocity Potential')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 平行一様流れの定義
def uniform_flow(x, y, U_inf):
return U_inf * np.ones_like(x), np.zeros_like(y)
# グリッドの作成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 流れの計算
U_inf = 1
U, V = uniform_flow(X, Y, U_inf)
# プロット
plt.figure(figsize=(10, 6))
plt.streamplot(X, Y, U, V, color='b')
plt.title('Uniform Flow')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# わき出しと吸い込みの定義
def source_sink_flow(x, y, strength):
r = np.sqrt(x**2 + y**2)
U = strength * x / (2 * np.pi * r**2)
V = -strength * y / (2 * np.pi * r**2)
return U, V
# グリッドの作成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 流れの計算
strength = 1
U, V = source_sink_flow(X, Y, strength)
# プロット
plt.figure(figsize=(10, 6))
plt.streamplot(X, Y, U, V, color='r')
plt.title('Source and Sink Flow')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 過流の定義
def doublet_flow(x, y, strength):
r2 = x**2 + y**2
U = -strength * x / (2 * np.pi * r2)
V = strength * y / (2 * np.pi * r2)
return U, V
# グリッドの作成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 流れの計算
strength = 1
U, V = doublet_flow(X, Y, strength)
# プロット
plt.figure(figsize=(10, 6))
plt.streamplot(X, Y, U, V, color='g')
plt.title('Doublet Flow')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 過流の定義
def doublet_flow(x, y, strength):
r2 = x**2 + y**2
U = -strength * x / (2 * np.pi * r2)
V = strength * y / (2 * np.pi * r2)
return U, V
# グリッドの作成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 流れの計算
strength = 1
U, V = doublet_flow(X, Y, strength)
# プロット
plt.figure(figsize=(10, 6))
plt.streamplot(X, Y, U, V, color='g')
plt.title('Doublet Flow')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 過流の定義
def doublet_flow(x, y, strength):
r2 = x**2 + y**2
U = -strength * x / (2 * np.pi * r2)
V = strength * y / (2 * np.pi * r2)
return U, V
# グリッドの作成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 流れの計算
strength = 1
U, V = doublet_flow(X, Y, strength)
# プロット
plt.figure(figsize=(10, 6))
plt.streamplot(X, Y, U, V, color='g')
plt.title('Doublet Flow')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 二重わき出しの定義
def doublet_source_flow(x, y, strength, a):
r2 = x**2 + y**2 + a**2
U = strength * x / (2 * np.pi * r2)
V = strength * y / (2 * np.pi * r2)
return U, V
# グリッドの作成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 流れの計算
strength = 1
a = 1
U, V = doublet_source_flow(X, Y, strength, a)
# プロット
plt.figure(figsize=(10, 6))
plt.streamplot(X, Y, U, V, color='b')
plt.title('Doublet with Source Flow')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 円柱まわりの流れの定義
def cylinder_flow(x, y, U_inf, R):
r2 = x**2 + y**2
U = U_inf * (1 - R**2 / r2)
V = -U_inf * (y / r2) * (R**2)
return U, V
# グリッドの作成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 円柱半径と流れの計算
R = 1
U_inf = 1
U, V = cylinder_flow(X, Y, U_inf, R)
# プロット
plt.figure(figsize=(10, 6))
plt.streamplot(X, Y, U, V, color='c')
plt.title('Flow Around a Cylinder')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
def boundary_layer_profile(x, delta_star, u_inf):
return u_inf * (1 - np.exp(-x / delta_star))
# パラメータ
x = np.linspace(0, 1, 100)
u_inf = 1
delta_star = 0.1
# 計算
u = boundary_layer_profile(x, delta_star, u_inf)
# プロット
plt.figure(figsize=(10, 6))
plt.plot(x, u, label='Boundary Layer Profile')
plt.title('Boundary Layer Profile')
plt.xlabel('x')
plt.ylabel('u')
plt.grid(True)
plt.show()
import numpy as np
# 抗力と揚力の計算
def forces(density, velocity, area, cd, cl):
drag = 0.5 * density * velocity**2 * area * cd
lift = 0.5 * density * velocity**2 * area * cl
return drag, lift
# パラメータ
density = 1.225 # kg/m^3
velocity = 10 # m/s
area = 1 # m^2
cd = 0.47 # 抗力係数
cl = 0.5 # 揚力係数
drag, lift = forces(density, velocity, area, cd, cl)
print(f'Drag Force: {drag} N')
print(f'Lift Force: {lift} N')
import numpy as np
import matplotlib.pyplot as plt
def karman_vortex(x, y, gamma):
r2 = x**2 + y**2
theta = np.arctan2(y, x)
u = -gamma / (2 * np.pi) * (np.sin(theta) / r2)
v = gamma / (2 * np.pi) * (np.cos(theta) / r2)
return u, v
# グリッドの作成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 流れの計算
gamma = 1.0 # 渦の強さ
U, V = karman_vortex(X, Y, gamma)
# プロット
plt.figure(figsize=(10, 6))
plt.streamplot(X, Y, U, V, color='r')
plt.title('Karman Vortex Street')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定義域
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 速度ポテンシャルと流れ関数の定義（均等流れの例）
U = 1 # 流れの速度
phi = U * X # 速度ポテンシャル
psi = U * Y # 流れ関数
# ベクトル場の計算
Vx, Vy = np.gradient(phi)
Vx, Vy = -Vy, Vx # 流れ関数から速度場を計算
# プロット
plt.figure(figsize=(10, 5))
# 速度ポテンシャルの等高線
plt.subplot(1, 2, 1)
plt.contourf(X, Y, phi, levels=20, cmap='viridis')
plt.colorbar(label='Velocity Potential')
plt.title('Velocity Potential')
# 流れ関数の等高線とベクトル場
plt.subplot(1, 2, 2)
plt.contour(X, Y, psi, levels=20, cmap='viridis')
plt.streamplot(X, Y, Vx, Vy, color='r', linewidth=1)
plt.title('Stream Function and Velocity Field')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定義域
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 複素ポテンシャルの定義（均等流れとソース流れの例）
U = 1 # 流れの速度
x0, y0 = 0, 0 # ソースの位置
W = U * X + 1j * U * Y + 1j * (np.log(np.sqrt((X - x0)**2 + (Y - y0)**2))) # 複素ポテンシャル
# 複素ポテンシャルの実部と虚部
phi = np.real(W)
psi = np.imag(W)
# ベクトル場の計算
Vx, Vy = np.gradient(phi)
Vx, Vy = -Vy, Vx # 流れ関数から速度場を計算
# プロット
plt.figure(figsize=(10, 5))
# 複素ポテンシャルの等高線
plt.subplot(1, 2, 1)
plt.contourf(X, Y, phi, levels=20, cmap='viridis')
plt.colorbar(label='Real Part of Complex Potential')
plt.title('Real Part of Complex Potential')
# 流れ関数の等高線とベクトル場
plt.subplot(1, 2, 2)
plt.contour(X, Y, psi, levels=20, cmap='viridis')
plt.streamplot(X, Y, Vx, Vy, color='r', linewidth=1)
plt.title('Stream Function and Velocity Field')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定義域
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 均等流れのポテンシャルと流れ関数
U = 1 # 流れの速度
phi = U * X
psi = U * Y
# プロット
plt.figure(figsize=(8, 4))
# 速度ポテンシャルの等高線
plt.subplot(1, 2, 1)
plt.contourf(X, Y, phi, levels=20, cmap='viridis')
plt.colorbar(label='Velocity Potential')
plt.title('Velocity Potential')
# 流れ関数の等高線
plt.subplot(1, 2, 2)
plt.contour(X, Y, psi, levels=20, cmap='viridis')
plt.title('Stream Function')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定義域
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 均等流れとソース流れのポテンシャル
U = 1 # 流れの速度
x0, y0 = 0, 0 # ソースの位置
phi = U * X + (1 / (2 * np.pi)) * np.log(np.sqrt((X - x0)**2 + (Y - y0)**2))
psi = U * Y - (1 / (2 * np.pi)) * np.arctan2(Y - y0, X - x0)
# プロット
plt.figure(figsize=(10, 5))
# 速度ポテンシャルの等高線
plt.subplot(1, 2, 1)
plt.contourf(X, Y, phi, levels=20, cmap='viridis')
plt.colorbar(label='Velocity Potential')
plt.title('Velocity Potential')
# 流れ関数の等高線
plt.subplot(1, 2, 2)
plt.contour(X, Y, psi, levels=20, cmap='viridis')
plt.title('Stream Function')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定義域
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 流れ場の定義（均等流れの例）
U = 1
phi = U * X # 速度ポテンシャル
# 圧力場の計算（簡略化された例）
pressure = -np.gradient(np.gradient(phi, axis=0), axis=0) - np.gradient(np.gradient(phi, axis=1), axis=1)
# プロット
plt.figure(figsize=(8, 4))
# 圧力の等高線
plt.contourf(X, Y, pressure, levels=20, cmap='viridis')
plt.colorbar(label='Pressure')
plt.title('Pressure Field')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# グリッドの生成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 複素数座標
Z = X + 1j * Y
# 簡単な流れの例: 平行流れ + 渦流れ
U = 1 - 0.5 * Y # 平行流れの成分
V = 0.5 * X # 渦流れの成分
# 流れ関数
psi = -0.5 * X**2 + 0.25 * Y**2
# 速度ポテンシャル
phi = X
# 複素速度ポテンシャル
w = phi + 1j * psi
# プロット
plt.figure(figsize=(12, 6))
# 流れ関数の等高線
plt.subplot(1, 2, 1)
plt.contour(X, Y, psi, levels=20, cmap='jet')
plt.title('Flow Function (Stream Function)')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar()
# 複素速度ポテンシャルの等高線
plt.subplot(1, 2, 2)
plt.contour(X, Y, np.real(w), levels=20, cmap='jet')
plt.title('Complex Velocity Potential (Real Part)')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar()
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
# ブラジウスの公式に基づく境界層の速度分布を求める関数
def blasius(f, y):
f1, f2, f3 = f
return [f2, f3, -0.5 * f1 * f3]
# 初期条件
f0 = [0, 0, 0] # [f(0), f'(0), f''(0)]
y_values = np.linspace(0, 10, 100) # yの範囲
# 数値的にブラジウスの方程式を解く
solution = odeint(blasius, f0, y_values)
f1 = solution[:, 0]
# 水面の微小振幅波の速度
g = 9.81 # 重力加速度 (m/s^2)
k = 2 * np.pi / 10 # 波数 (1/m)
c = np.sqrt(g / k) # 波速度 (m/s)
# プロット
plt.figure(figsize=(12, 6))
# ブラジウスの境界層速度分布
plt.subplot(1, 2, 1)
plt.plot(y_values, f1, label='Boundary Layer Velocity Profile')
plt.xlabel('η')
plt.ylabel('f(η)')
plt.title('Blasius Boundary Layer Velocity Profile')
plt.grid(True)
plt.legend()
# 水面の波速度
plt.subplot(1, 2, 2)
plt.bar(['Wave Speed'], [c], color='skyblue')
plt.ylabel('Speed (m/s)')
plt.title('Wave Speed for Small Amplitude Waves')
plt.grid(True)
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# グリッドの生成
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)
# 速度場の定義（単純なせん断流れの例）
U = -Y
V = X
# 速度場の勾配の計算
du_dx, du_dy = np.gradient(U, x, y)
dv_dx, dv_dy = np.gradient(V, x, y)
# 変形速度テンソルの計算
epsilon_xx = du_dx
epsilon_yy = dv_dy
epsilon_xy = 0.5 * (du_dy + dv_dx)
# 応力テンソルの計算（動粘性係数を仮定）
mu = 1.0 # 動粘性係数
tau_xx = mu * epsilon_xx
tau_yy = mu * epsilon_yy
tau_xy = mu * epsilon_xy
# プロット
plt.figure(figsize=(12, 6))
# 変形速度テンソルのプロット
plt.subplot(1, 2, 1)
plt.contourf(X, Y, epsilon_xy, cmap='jet')
plt.colorbar(label='Shear Strain Rate (s$^{-1}$)')
plt.title('Shear Strain Rate')
plt.xlabel('x')
plt.ylabel('y')
# 応力テンソルのプロット
plt.subplot(1, 2, 2)
plt.contourf(X, Y, tau_xy, cmap='jet')
plt.colorbar(label='Shear Stress (Pa)')
plt.title('Shear Stress')
plt.xlabel('x')
plt.ylabel('y')
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# ベクトルの演算
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
# ベクトルの加算
sum_ab = a + b
# ベクトルのスカラー倍
c = 2
scaled_a = c * a
# ベクトルの内積
dot_product = np.dot(a, b)
# ベクトルの外積
cross_product = np.cross(a, b)
# ベクトル場の微分
# スカラー場とベクトル場の例
def gradient(f, x, y):
fx, fy = np.gradient(f, x, y)
return fx, fy
# スカラー場の定義
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)
Z = X**2 + Y**2 # スカラー場
# 勾配の計算
fx, fy = gradient(Z, x, y)
# ベクトル場の積分
# 簡単な線積分の例
def line_integral(vector_field, path):
integral = 0
for i in range(len(path) - 1):
dx = path[i + 1] - path[i]
integral += np.dot(vector_field[i], dx)
return integral
# 簡単なベクトル場と経路
vector_field = np.array([[1, 0] for _ in range(len(x))])
path = np.array(x)
# 線積分の計算
line_integral_result = line_integral(vector_field, path)
# テンソルの例
tensor = np.array([[1, 2], [3, 4]])
scalar = 2
scaled_tensor = scalar * tensor
# プロット
plt.figure(figsize=(12, 6))
# スカラー場の勾配
plt.subplot(1, 2, 1)
plt.contourf(X, Y, Z, cmap='viridis')
plt.quiver(X, Y, fx, fy, color='red')
plt.title('Scalar Field and Gradient')
plt.colorbar()
# テンソルのスカラー倍
plt.subplot(1, 2, 2)
plt.imshow(scaled_tensor, cmap='viridis', interpolation='none')
plt.title('Scaled Tensor')
plt.colorbar()
plt.tight_layout()
plt.show()
R = 1.0 # 円柱の半径
# 円柱まわりの流れの計算
W_cylinder = U * (Z + R**2 / Z)
# 流線の計算
streamlines_cylinder = np.real(W_cylinder)
contour_lines_cylinder = np.imag(W_cylinder)
# 可視化
plt.figure(figsize=(10, 6))
plt.contour(X, Y, streamlines_cylinder, levels=20, colors='blue')
plt.contour(X, Y, contour_lines_cylinder, levels=20, colors='red')
plt.title('Streamlines around a Cylinder')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
alpha = 1.0 # わき出しまたは吸込みの強さ
# わき出しと吸込みの影響の計算
W_sink_source = U * (Z - alpha / Z)
# 流線の計算
streamlines_sink_source = np.real(W_sink_source)
contour_lines_sink_source = np.imag(W_sink_source)
# 可視化
plt.figure(figsize=(10, 6))
plt.contour(X, Y, streamlines_sink_source, levels=20, colors='blue')
plt.contour(X, Y, contour_lines_sink_source, levels=20, colors='red')
plt.title('Streamlines with Sink and Source')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
# 複素ポテンシャルの計算
Z = X + 1j * Y
W = U * Z
# 流線の計算
streamlines = np.real(W)
contour_lines = np.imag(W)
# 可視化
plt.figure(figsize=(10, 6))
plt.contour(X, Y, streamlines, levels=20, colors='blue')
plt.contour(X, Y, contour_lines, levels=20, colors='red')
plt.title('Streamlines of Uniform Flow using Complex Potential')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
# 流れ関数の計算
Psi = U * Y
# 可視化
plt.figure(figsize=(10, 6))
plt.contourf(X, Y, Psi, cmap='coolwarm', levels=100)
plt.colorbar(label='Stream Function (Ψ)')
plt.title('Uniform Flow Stream Function')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定数の設定
U = 1.0 # 流れの速度
x = np.linspace(-10, 10, 100)
y = np.linspace(-10, 10, 100)
X, Y = np.meshgrid(x, y)
# 速度ポテンシャルの計算
Phi = U * X
# 可視化
plt.figure(figsize=(10, 6))
plt.contourf(X, Y, Phi, cmap='coolwarm', levels=100)
plt.colorbar(label='Velocity Potential (Φ)')
plt.title('Uniform Flow Velocity Potential')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
Gamma = 1.0 # 渦度
# 角を回る流れの計算
W_vortex = Gamma * np.log(Z)
# 流線の計算
streamlines_vortex = np.real(W_vortex)
contour_lines_vortex = np.imag(W_vortex)
# 可視化
plt.figure(figsize=(10, 6))
plt.contour(X, Y, streamlines_vortex, levels=20, colors='blue')
plt.contour(X, Y, contour_lines_vortex, levels=20, colors='red')
plt.title('Streamlines of Vortex Flow')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
alpha1 = 1.0
alpha2 = 0.5
alpha3 = -0.5
z1 = 2.0
z2 = -2.0
# 多重のわき出し吸い込みの計算
W_multiple_sources = U * (Z + alpha1 / Z + alpha2 / (Z - z1) + alpha3 / (Z - z2))
# 流線の計算
streamlines_multiple_sources = np.real(W_multiple_sources)
contour_lines_multiple_sources = np.imag(W_multiple_sources)
# 可視化
plt.figure(figsize=(10, 6))
plt.contour(X, Y, streamlines_multiple_sources, levels=20, colors='blue')
plt.contour(X, Y, contour_lines_multiple_sources, levels=20, colors='red')
plt.title('Streamlines of Multiple Sources/Sinks')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
# ジューコフスキー変換
R = 1.0
zeta = Z + R**2 / Z
# 流線の計算
streamlines_zeta = np.real(zeta)
contour_lines_zeta = np.imag(zeta)
# 可視化
plt.figure(figsize=(10, 6))
plt.contour(X, Y, streamlines_zeta, levels=20, colors='blue')
plt.contour(X, Y, contour_lines_zeta, levels=20, colors='red')
plt.title('Streamlines after Joukowski Transformation')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 流れの速度場を定義
def velocity_field(x, y):
u = -y
v = x
return u, v
# 流量の計算（簡略化されたモデル）
def flow_rate(u, v, area):
return np.sum(u * area[0] + v * area[1])
x = np.linspace(-1, 1, 10)
y = np.linspace(-1, 1, 10)
X, Y = np.meshgrid(x, y)
U, V = velocity_field(X, Y)
# 流量を計算
area = np.ones_like(U) # 単位断面積
Q = flow_rate(U, V, area)
print(f'流量 Q = {Q}')
# 加速度の計算（簡略化されたモデル）
def acceleration(u, v, dx, dy):
du_dx = np.gradient(u, axis=1) / dx
du_dy = np.gradient(u, axis=0) / dy
dv_dx = np.gradient(v, axis=1) / dx
dv_dy = np.gradient(v, axis=0) / dy
a_x = du_dx + u * du_dx + v * dv_dx
a_y = du_dy + u * du_dy + v * dv_dy
return a_x, a_y
dx, dy = x[1] - x[0], y[1] - y[0]
A_x, A_y = acceleration(U, V, dx, dy)
# 可視化
plt.figure(figsize=(10, 6))
plt.quiver(X, Y, A_x, A_y, scale=10, color='r')
plt.title('Acceleration Field')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
# 流線のプロット
plt.figure(figsize=(10, 6))
plt.streamplot(X, Y, U, V, color='b', linewidth=2)
plt.title('Streamlines')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
# 定常流れの例
def steady_flow(x, y):
return -y, x
# 非定常流れの例
def unsteady_flow(x, y, t):
return -y * np.sin(t), x * np.cos(t)
# 時間ステップ
t = 1.0
# 流れの計算
U_steady, V_steady = steady_flow(X, Y)
U_unsteady, V_unsteady = unsteady_flow(X, Y, t)
# 可視化
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.streamplot(X, Y, U_steady, V_steady, color='b', linewidth=2)
plt.title('Steady Flow')
plt.subplot(1, 2, 2)
plt.streamplot(X, Y, U_unsteady, V_unsteady, color='r', linewidth=2)
plt.title('Unsteady Flow')
plt.show()
# 一様流れ
U_uniform = np.ones_like(X)
V_uniform = np.zeros_like(Y)
# 非一様流れ
def nonuniform_flow(x, y):
return -y * np.exp(-x), x * np.exp(-y)
U_nonuniform, V_nonuniform = nonuniform_flow(X, Y)
# 可視化
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.streamplot(X, Y, U_uniform, V_uniform, color='b', linewidth=2)
plt.title('Uniform Flow')
plt.subplot(1, 2, 2)
plt.streamplot(X, Y, U_nonuniform, V_nonuniform, color='r', linewidth=2)
plt.title('Nonuniform Flow')
plt.show()
Gamma = 1.0
def vortex_flow(x, y):
r2 = x**2 + y**2
u = -Gamma / (2 * np.pi) * (y / r2)
v = Gamma / (2 * np.pi) * (x / r2)
return u, v
U_vortex, V_vortex = vortex_flow(X, Y)
# 可視化
plt.figure(figsize=(10, 6))
plt.streamplot(X, Y, U_vortex, V_vortex, color='purple', linewidth=2)
plt.title('Vortex Flow')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
# 混相流の簡単なモデル
def multiphase_flow(x, y, t):
return -y * np.sin(t), x * np.cos(t)
U_multiphase, V_multiphase = multiphase_flow(X, Y, t)
# 可視化
plt.figure(figsize=(10, 6))
plt.streamplot(X, Y, U_multiphase, V_multiphase, color='orange', linewidth=2)
plt.title('Multiphase Flow')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 速度場の定義
def velocity_field(x, y):
# 例: 一様流れ
v_x = np.ones_like(x)
v_y = np.zeros_like(y)
return v_x, v_y
# グリッドの生成
x = np.linspace(0, 10, 20)
y = np.linspace(0, 10, 20)
X, Y = np.meshgrid(x, y)
V_x, V_y = velocity_field(X, Y)
# 速度場のプロット
plt.figure(figsize=(8, 6))
plt.quiver(X, Y, V_x, V_y, scale=5)
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Velocity Field')
plt.grid()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
# グリッドの設定
nx, ny = 100, 100
x = np.linspace(0, 2*np.pi, nx)
y = np.linspace(0, 2*np.pi, ny)
X, Y = np.meshgrid(x, y)
# 速度場の定義 (非圧縮性流体の例)
# 速度場の定義を変更して壁や障害物を考慮することができます
U = -np.sin(X) * np.cos(Y)
V = np.cos(X) * np.sin(Y)
# ダイバージェンスの計算
div_U = np.gradient(U, x, axis=0)
div_V = np.gradient(V, y, axis=1)
div = div_U + div_V
# 可視化
plt.figure(figsize=(14, 7))
# 速度ベクトル場のプロット
plt.subplot(1, 2, 1)
plt.quiver(X, Y, U, V, scale=5, color='b')
plt.title('Velocity Field')
plt.xlabel('x')
plt.ylabel('y')
plt.xlim([0, 2*np.pi])
plt.ylim([0, 2*np.pi])
# ダイバージェンスのプロット
plt.subplot(1, 2, 2)
plt.contourf(X, Y, div, levels=20, cmap='RdBu_r', norm=Normalize(vmin=-0.4, vmax=0.4))
plt.colorbar(label='Divergence')
plt.title('Divergence of Velocity Field')
plt.xlabel('x')
plt.ylabel('y')
plt.xlim([0, 2*np.pi])
plt.ylim([0, 2*np.pi])
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定数の定義
mu = 0.001 # 粘度 (Pa.s)
v = 0.1 # 上部ガラス板の速度 (m/s)
d = 0.01 # 流体の厚み (m)
# グリッドの設定
y = np.linspace(0, d, 100)
v_x = np.zeros_like(y)
# 速度の計算 (クエット流れ)
v_x = (v / d) * y
# 接線応力の計算
tau = mu * (v / d) * np.ones_like(y)
# プロット設定
fig, ax = plt.subplots(2, 1, figsize=(10, 8))
# 速度プロファイルのプロット
ax[0].plot(v_x, y, label='Velocity Profile')
ax[0].set_xlabel('Velocity (m/s)')
ax[0].set_ylabel('Distance from bottom plate (m)')
ax[0].set_title('Velocity Profile in Couette Flow')
ax[0].grid(True)
ax[0].legend()
# 接線応力プロファイルのプロット
ax[1].plot(tau, y, label='Shear Stress', color='r')
ax[1].set_xlabel('Shear Stress (Pa)')
ax[1].set_ylabel('Distance from bottom plate (m)')
ax[1].set_title('Shear Stress in Couette Flow')
ax[1].grid(True)
ax[1].legend()
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import laplace
# シミュレーションのパラメータ
nx, ny = 100, 100 # グリッドサイズ
dx = dy = 1.0 # グリッド間隔
dt = 0.01 # 時間ステップ
nu = 0.1 # 動粘性率 μ
lambda_ = 0.2 # 体積粘性率 λ
iterations = 100 # シミュレーションの反復回数
# 初期速度場
u = np.zeros((nx, ny))
v = np.zeros((nx, ny))
# 初期密度場 (中心に高密度)
rho = np.ones((nx, ny))
rho[nx//2-5:nx//2+5, ny//2-5:ny//2+5] = 2.0
def divergence(u, v, dx, dy):
"""速度場の発散を計算"""
du_dx = (np.roll(u, -1, axis=0) - np.roll(u, 1, axis=0)) / (2 * dx)
dv_dy = (np.roll(v, -1, axis=1) - np.roll(v, 1, axis=1)) / (2 * dy)
return du_dx + dv_dy
def update_velocity(u, v, rho, nu, lambda_, dx, dy, dt):
"""速度場を更新"""
div_v = divergence(u, v, dx, dy)
# 圧力勾配（簡略化）
pressure_grad_u = -np.gradient(rho, axis=0)
pressure_grad_v = -np.gradient(rho, axis=1)
# ラプラシアンによる粘性項
lap_u = laplace(u)
lap_v = laplace(v)
# 速度場の更新
u += dt * (-u * np.gradient(u, axis=0) - v * np.gradient(u, axis=1)
+ nu * lap_u + lambda_ * div_v + pressure_grad_u)
v += dt * (-u * np.gradient(v, axis=0) - v * np.gradient(v, axis=1)
+ nu * lap_v + lambda_ * div_v + pressure_grad_v)
# シミュレーション
for _ in range(iterations):
update_velocity(u, v, rho, nu, lambda_, dx, dy, dt)
# 結果の可視化
plt.quiver(u, v)
plt.title('Velocity Field')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import laplace
# シミュレーションのパラメータ
nx, ny = 100, 100 # グリッドサイズ
dx = dy = 1.0 # グリッド間隔
dt = 0.01 # 時間ステップ
nu = 0.1 # 動粘性率 μ
alpha = 0.01 # 熱拡散率
gamma = 1.4 # 比熱比
C = 1.0 # 状態方程式の定数
iterations = 500 # シミュレーションの反復回数
# 初期条件設定
u = np.zeros((nx, ny))
v = np.zeros((nx, ny))
T = np.ones((nx, ny)) * 300 # 温度を300Kに初期化
T[nx//2-5:nx//2+5, ny//2-5:ny//2+5] = 600 # 中央に局所的な高温領域を設定
def divergence(u, v, dx, dy):
"""速度場の発散を計算"""
du_dx = (np.roll(u, -1, axis=0) - np.roll(u, 1, axis=0)) / (2 * dx)
dv_dy = (np.roll(v, -1, axis=1) - np.roll(v, 1, axis=1)) / (2 * dy)
return du_dx + dv_dy
def update_temperature(T, u, v, alpha, dx, dy, dt):
"""温度場を更新"""
T_new = T.copy()
T_new += alpha * laplace(T) * dt # 熱拡散
T_new -= u * (np.roll(T, -1, axis=0) - np.roll(T, 1, axis=0)) / (2 * dx) * dt # 温度の対流項
T_new -= v * (np.roll(T, -1, axis=1) - np.roll(T, 1, axis=1)) / (2 * dy) * dt
return T_new
def update_velocity(u, v, T, nu, gamma, C, dx, dy, dt):
"""速度場を更新"""
rho = C / T # 状態方程式の仮定: 密度は温度に反比例
p = rho * T # 理想気体の状態方程式から圧力を計算
div_v = divergence(u, v, dx, dy)
# 圧力勾配
pressure_grad_u = -np.gradient(p, axis=0)
pressure_grad_v = -np.gradient(p, axis=1)
# ラプラシアンによる粘性項
lap_u = laplace(u)
lap_v = laplace(v)
# 速度場の更新
u += dt * (-u * np.gradient(u, axis=0) - v * np.gradient(u, axis=1)
+ nu * lap_u - pressure_grad_u / rho)
v += dt * (-u * np.gradient(v, axis=0) - v * np.gradient(v, axis=1)
+ nu * lap_v - pressure_grad_v / rho)
return u, v
# シミュレーションのメインループ
for _ in range(iterations):
T = update_temperature(T, u, v, alpha, dx, dy, dt)
u, v = update_velocity(u, v, T, nu, gamma, C, dx, dy, dt)
# 結果の可視化
plt.figure(figsize=(12, 5))
# 温度分布のプロット
plt.subplot(1, 2, 1)
plt.imshow(T, cmap='hot', origin='lower')
plt.colorbar()
plt.title('Temperature Field')
# 速度場のプロット
plt.subplot(1, 2, 2)
plt.quiver(u, v)
plt.title('Velocity Field')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import laplace
# グリッドサイズとパラメータ
nx, ny = 100, 100 # グリッドサイズ
dx, dy = 1.0, 1.0 # グリッド間隔
dt = 0.01 # 時間ステップ
nu = 0.1 # 動粘性率 μ
rho = 1.0 # 密度 ρ
iterations = 100 # シミュレーションの反復回数
# 初期条件の設定
u = np.zeros((nx, ny)) # x方向の速度場
v = np.zeros((nx, ny)) # y方向の速度場
p = np.zeros((nx, ny)) # 圧力場
# 圧力差を与えるための初期設定
p[:, ny//2:] = 1.0 # 上部に高圧、下部に低圧を設定
def update_velocity(u, v, p, rho, nu, dx, dy, dt):
"""速度場の更新"""
lap_u = laplace(u) / (dx * dy)
lap_v = laplace(v) / (dx * dy)
# 圧力勾配
dpdx = np.gradient(p, axis=1) / dx
dpdy = np.gradient(p, axis=0) / dy
# Navier-Stokes方程式の簡略化
u_new = u + dt * (-u * np.gradient(u, axis=1) - v * np.gradient(u, axis=0)
+ nu * lap_u - dpdx / rho)
v_new = v + dt * (-u * np.gradient(v, axis=1) - v * np.gradient(v, axis=0)
+ nu * lap_v - dpdy / rho)
return u_new, v_new
def update_pressure(p, u, v, dx, dy, dt):
"""圧力場の更新"""
div_v = np.gradient(u, axis=1) / dx + np.gradient(v, axis=0) / dy
p_new = p - dt * rho * div_v # 簡略化された圧力更新
return p_new
# メインのシミュレーションループ
for _ in range(iterations):
u, v = update_velocity(u, v, p, rho, nu, dx, dy, dt)
p = update_pressure(p, u, v, dx, dy, dt)
# 結果の可視化
plt.figure(figsize=(12, 5))
# 圧力場のプロット
plt.subplot(1, 2, 1)
plt.imshow(p, cmap='jet', origin='lower')
plt.colorbar()
plt.title('Pressure Field')
# 速度場のプロット
plt.subplot(1, 2, 2)
plt.quiver(u, v)
plt.title('Velocity Field')
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# 定数
g = 9.81 # 重力加速度 (m/s^2)
rho = 1.0 # 流体の密度 (kg/m^3)
# 流体の状態を定義
# ここでは単位体積あたりのエネルギーを計算するためのサンプルデータを作成します
heights = np.linspace(0, 10, 100) # 高さ (m)
velocities = np.linspace(0, 20, 100) # 速度 (m/s)
pressures = np.linspace(100000, 1000000, 100) # 圧力 (Pa)
# ベルヌーイの定理を計算
def bernoulli_theorem(v, p, h):
return 0.5 * v**2 + p / rho + g * h
# 高さに対するエネルギーをプロットする
def plot_energy():
plt.figure(figsize=(12, 6))
# 速度と圧力が固定の場合のエネルギー
v_fixed = 10 # 固定速度 (m/s)
p_fixed = 500000 # 固定圧力 (Pa)
energy = bernoulli_theorem(v_fixed, p_fixed, heights)
plt.plot(heights, energy, label=f'Fixed velocity {v_fixed} m/s, pressure {p_fixed} Pa')
# 高さに対するエネルギーをプロット
for v, p in zip(velocities, pressures):
energy = bernoulli_theorem(v, p, heights)
plt.plot(heights, energy, label=f'v = {v:.1f} m/s, p = {p:.0f} Pa')
plt.xlabel('Height (m)')
plt.ylabel('Total Energy per Unit Mass (J/kg)')
plt.title('Bernoulli\'s Theorem: Energy vs. Height')
plt.legend()
plt.grid()
plt.show()
# プロットを作成
plot_energy()
import numpy as np
import matplotlib.pyplot as plt
# パラメータ
Gamma = 1.0 # 渦強度 (例)
U_inf = 1.0 # 気流の速度 (例)
r_max = 2.0 # 半径の最大値
resolution = 100 # 解像度
# 半径と角度のグリッドを生成
r = np.linspace(0.01, r_max, resolution) # r=0を避けるため、0.01から開始
theta = np.linspace(0, 2 * np.pi, resolution)
R, Theta = np.meshgrid(r, theta)
# 渦の速度場の計算
U_vortex = Gamma / (2 * np.pi * R)
U_vortex_x = U_vortex * np.cos(Theta)
U_vortex_y = U_vortex * np.sin(Theta)
# 気流の速度場の計算
U_flow_x = U_inf * np.ones_like(R)
U_flow_y = np.zeros_like(R)
# 渦と気流の速度場を合成
U_total_x = U_flow_x + U_vortex_x
U_total_y = U_flow_y + U_vortex_y
# Cartesian座標への変換
X = R * np.cos(Theta)
Y = R * np.sin(Theta)
# 可視化
plt.figure(figsize=(12, 6))
# 渦の速度場
plt.subplot(1, 2, 1)
plt.quiver(X, Y, U_vortex_x, U_vortex_y, scale=5, color='r')
plt.title('Vortex Flow Field')
plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.axis('equal')
plt.grid()
# 気流と渦の合成速度場
plt.subplot(1, 2, 2)
plt.quiver(X, Y, U_total_x, U_total_y, scale=5, color='b')
plt.title('Combined Flow Field')
plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.axis('equal')
plt.grid()
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
# レイノルズ数の計算
def reynolds_number(density, velocity, diameter, viscosity):
return (density * velocity * diameter) / viscosity
# 圧力損失の計算 (層流用)
def pressure_loss_laminar(flow_rate, length, viscosity, diameter):
return (8 * viscosity * flow_rate * length) / (np.pi * diameter**4)
# 圧力損失の計算 (乱流用)
def pressure_loss_turbulent(flow_rate, length, diameter, density, friction_factor):
velocity = flow_rate / (np.pi * (diameter/2)**2)
return friction_factor * (length/diameter) * (density * velocity**2 / 2)
# ダルシー・ワイスバッハ摩擦係数の計算
def darcy_weisbach_friction_factor(Re, roughness, diameter):
if Re < 2300:
return 64 / Re
else:
def colebrook(f):
return 1/np.sqrt(f) + 2.0*np.log10(roughness/(3.7*diameter) + 2.51/(Re*np.sqrt(f)))
initial_guess = 0.02
return fsolve(colebrook, initial_guess)[0]
# 管路抵抗の計算
def pipe_resistance(pressure_loss, flow_rate):
return pressure_loss / flow_rate
# パラメータ設定
density = 1000 # kg/m^3
velocity = 2.0 # m/s
diameter = 0.05 # m
viscosity = 0.001 # Pa.s
flow_rate = 0.01 # m^3/s
length = 10.0 # m
roughness = 0.0001 # m
# 各種計算
Re = reynolds_number(density, velocity, diameter, viscosity)
pressure_loss = pressure_loss_laminar(flow_rate, length, viscosity, diameter)
friction_factor = darcy_weisbach_friction_factor(Re, roughness, diameter)
pressure_loss_turb = pressure_loss_turbulent(flow_rate, length, diameter, density, friction_factor)
pipe_resistance_value = pipe_resistance(pressure_loss_turb, flow_rate)
# 結果出力
print(f"レイノルズ数: {Re}")
print(f"層流の圧力損失: {pressure_loss} Pa")
print(f"ダルシー・ワイスバッハ摩擦係数: {friction_factor}")
print(f"乱流の圧力損失: {pressure_loss_turb} Pa")
print(f"管路抵抗: {pipe_resistance_value} Pa.s/m^3")