流体力学Pythonコード



```
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")
```


```


```


```