模拟山体滑坡

最新推荐文章于 2025-05-01 23:36:42 发布

骑蜗牛上月亮

最新推荐文章于 2025-05-01 23:36:42 发布

阅读量1.2k

点赞数 27

文章标签： python

本文链接：https://blog.csdn.net/qq_42788199/article/details/147620402

版权

模拟山体滑坡：利用简单的物理模型和 matplotlib，模拟山体上不同质量的物体在一定坡度下的滑动过程，展示滑坡的形成。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.ndimage import gaussian_filter


class LandslideSimulator:
    def __init__(self):
        self.grid_size = (200, 200)
        self.dx = 1.0
        self.dt = 0.1
        self.g = 9.81

        self.density = 1600.0
        self.cohesion = 5000.0
        self.friction_angle = 30.0
        self.pore_pressure = 0.3

        self.initialize_terrain()

        self.velocity = np.zeros(self.grid_size + (2,))
        self.stress = np.zeros(self.grid_size + (2, 2))

    def initialize_terrain(self):
        x = np.linspace(0, 20, self.grid_size[1])
        y = np.linspace(0, 20, self.grid_size[0])
        X, Y = np.meshgrid(x, y)
        self.height = 15 - 0.4 * X + 2 * np.sin(0.5 * X)

        noise = np.random.normal(0, 0.5, self.grid_size)
        self.height += gaussian_filter(noise, sigma=2)

        self.strength = np.ones(self.grid_size) * self.cohesion
        weakness_zone = (X > 8) & (X < 12) & (Y > 5) & (Y < 15)
        self.strength[weakness_zone] *= 0.2

        self.mask = np.ones(self.grid_size, dtype=bool)
        self.original_height = self.height.copy()

    def calculate_stress(self):
        depth = self.original_height - self.height
        # 检查深度值，避免数值过大
        depth = np.clip(depth, 0, 10)
        sigma_z = self.density * self.g * depth

        tau_max = (sigma_z * (1 - self.pore_pressure) *
                   np.tan(np.deg2rad(self.friction_angle)) + self.strength)

        grad = np.gradient(self.height, self.dx)
        tau_x = -self.density * self.g * grad[0]
        tau_y = -self.density * self.g * grad[1]

        self.stress[..., 0, 0] = tau_x
        self.stress[..., 1, 1] = tau_y
        self.stress[..., 0, 1] = self.stress[..., 1, 0] = 0.5 * (tau_x + tau_y)

        # 修正后的计算表达式（原71行）
        tau_magnitude = np.sqrt(tau_x ** 2 + tau_y ** 2)
        self.failure_mask = tau_magnitude > tau_max

    def update_velocity(self):
        body_force = np.zeros_like(self.velocity)
        body_force[..., 1] = self.density * self.g

        stress_grad = np.gradient(self.stress, self.dx, axis=(0, 1))
        total_force = body_force - (stress_grad[0][..., 0] + stress_grad[1][..., 1])

        # 检查总力值，避免数值过大
        total_force = np.clip(total_force, -1e6, 1e6)
        acceleration = total_force / self.density
        self.velocity += acceleration * self.dt
        self.velocity *= 0.98

    def update_height(self):
        flux = np.zeros_like(self.height)

        vx = self.velocity[..., 0]
        height_diff_x = self.height[:, :-2] - self.height[:, 2:]
        # 检查高度差，避免数值过大
        height_diff_x = np.clip(height_diff_x, -1e6, 1e6)
        flux[:, 1:-1] += vx[:, 1:-1] * height_diff_x / (2 * self.dx)

        vy = self.velocity[..., 1]
        height_diff_y = self.height[:-2, :] - self.height[2:, :]
        # 检查高度差，避免数值过大
        height_diff_y = np.clip(height_diff_y, -1e6, 1e6)
        flux[1:-1, :] += vy[1:-1, :] * height_diff_y / (2 * self.dx)

        # 检查通量值，避免数值过大
        flux = np.clip(flux, -1e6, 1e6)
        self.height += flux * self.dt
        self.height[self.failure_mask] -= 0.1 * self.dt * self.velocity[self.failure_mask, 1]

        self.height[:, 0] = self.original_height[:, 0]
        self.height[0, :] = self.original_height[0, :]

    def step(self):
        self.calculate_stress()
        self.update_velocity()
        self.update_height()


# 初始化模拟器和图形
sim = LandslideSimulator()
fig = plt.figure(figsize=(10, 5))
ax1 = fig.add_subplot(121, projection='3d')
ax2 = fig.add_subplot(122)

x = np.linspace(0, 20, sim.grid_size[1])
y = np.linspace(0, 20, sim.grid_size[0])
X, Y = np.meshgrid(x, y)


def update(frame):
    sim.step()
    ax1.cla()
    ax2.cla()

    ax1.plot_surface(X, Y, sim.height, cmap='terrain', alpha=0.8)
    ax1.set_zlim(0, 20)
    ax1.set_title('3D opographic Evolution')

    skip = 5
    ax2.quiver(X[::skip, ::skip], Y[::skip, ::skip],
               sim.velocity[::skip, ::skip, 0],
               sim.velocity[::skip, ::skip, 1],
               scale=50, color='red')
    ax2.imshow(sim.height, extent=[0, 20, 0, 20], cmap='terrain', origin='lower')

    return ax1, ax2


ani = FuncAnimation(fig, update, frames=200, interval=50, blit=False)
plt.show()