生成一个20*35的矩形 其为一容器侧视图 可绕其左下角顶点旋转 旋转角度为0°到180° 矩形装满水 水会随着矩形旋转而流出 通过控制 矩形的旋转角速度 使水流量保持恒定
import math
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
class EnhancedRotatingRectangle:
def __init__(self, width, height):
self.W = width
self.H = height
self.alpha = math.atan(height / width) # 临界角度(弧度)
self.full_area = width * height
# 添加角度限制参数
self.min_angle = 0
self.max_angle = math.pi
# 添加流体特性参数
self.flow_coefficient = 0.62 # 流量系数(考虑粘滞性)
def area(self, theta):
"""优化后的水量计算(包含90°后的状态)"""
if theta < self.min_angle:
return self.full_area
elif theta <= self.alpha:
return self.W * self.H - 0.5 * self.W**2 * math.tan(theta)
elif theta < math.pi/2:
return 0.5 * self.H**2 / math.tan(theta)
elif theta < math.pi:
# 倒置状态的水量计算
return 0.5 * self.W**2 * math.tan(math.pi - theta)
else:
return 0.0
def dA_dtheta(self, theta):
"""水量对旋转角的导数(高阶精度)"""
if theta < 0 or theta > math.pi:
return 0
# 使用中心差分法提高精度
epsilon = 1e-6
A_plus = self.area(theta + epsilon)
A_minus = self.area(theta - epsilon)
return (A_plus - A_minus) / (2 * epsilon)
def angular_velocity(self, theta, flow_rate):
"""改进的角速度计算(包含流体特性)"""
if theta < self.min_angle or theta > self.max_angle:
return 0.0
dA_dtheta = self.dA_dtheta(theta)
# 考虑流量系数
adjusted_flow_rate = flow_rate / self.flow_coefficient
if dA_dtheta >= 0:
return 0.0
return adjusted_flow_rate / (-dA_dtheta)
def simulate_rotation(self, flow_rate, dt=0.01):
"""封装模拟过程并返回结果"""
total_time = 0.0
theta_rad = 0.0
results = {
'time': [0.0],
'theta_rad': [0.0],
'omega': [self.angular_velocity(0, flow_rate)],
'area': [self.area(0)],
'flow_rate_actual': [0.0]
}
# 从0°模拟到180°
while theta_rad < self.max_angle:
# 计算角速度
omega = self.angular_velocity(theta_rad, flow_rate)
# 角度更新(避免超出范围)
d_theta = omega * dt
if theta_rad + d_theta > self.max_angle:
theta_rad = self.max_angle
else:
theta_rad += d_theta
total_time += dt
# 记录数据
results['time'].append(total_time)
results['theta_rad'].append(theta_rad)
results['omega'].append(omega)
results['area'].append(self.area(theta_rad))
# 计算实际流量
if len(results['time']) > 1:
dA = results['area'][-2] - results['area'][-1]
dt_step = results['time'][-1] - results['time'][-2]
flow_actual = dA / dt_step
else:
flow_actual = 0
results['flow_rate_actual'].append(flow_actual)
return results
# 参数设置(与原始代码一致)
W = 20.0 # 宽度
H = 35.0 # 高度
Q = 10.0 # 恒定流量
dt = 0.01 # 时间步长
# 创建增强版模型实例
enhanced_rect = EnhancedRotatingRectangle(W, H)
# 运行模拟
sim_results = enhanced_rect.simulate_rotation(Q, dt)
# 可视化结果
plt.figure(figsize=(14, 12))
# 1. 角度变化
plt.subplot(411)
plt.plot(sim_results['time'], [math.degrees(t) for t in sim_results['theta_rad']])
plt.ylabel('Rotation Angle (°)')
plt.title('Container Rotation and Water Flow Control')
plt.grid(True)
# 2. 角速度变化
plt.subplot(412)
plt.plot(sim_results['time'], [math.degrees(w) for w in sim_results['omega']])
plt.ylabel('Angular Velocity (°/s)')
plt.grid(True)
# 3. 水量变化
plt.subplot(413)
plt.plot(sim_results['time'], sim_results['area'])
plt.ylabel('Water Area (m²)')
plt.grid(True)
# 4. 流量对比
plt.subplot(414)
plt.plot(sim_results['time'], sim_results['flow_rate_actual'], label='Actual Flow')
plt.axhline(y=Q, color='r', linestyle='--', label='Target Flow')
plt.ylabel('Flow Rate (m²/s)')
plt.xlabel('Time (s)')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig('water_flow_control.png', dpi=300)
plt.show()
# 创建动态可视化
fig, ax = plt.subplots(figsize=(10, 8))
ax.set_xlim(-5, 40)
ax.set_ylim(-5, 40)
ax.set_aspect('equal')
ax.grid(True)
ax.set_title('Rotating Container with Water')
# 创建容器和水面对象
container_line, = ax.plot([], [], 'b-', linewidth=2)
water_poly, = ax.fill([], [], color='cyan', alpha=0.6)
pivot_point = ax.plot(0, 0, 'ro', markersize=8)[0]
theta_text = ax.text(0.02, 0.95, '', transform=ax.transAxes)
flow_text = ax.text(0.02, 0.90, '', transform=ax.transAxes)
def init():
container_line.set_data([], [])
water_poly.set_xy(np.empty((0, 2)))
return container_line, water_poly, pivot_point
def update(frame):
# 获取当前状态
theta = sim_results['theta_rad'][frame]
area = sim_results['area'][frame]
# 计算容器角点
corners = np.array([
[0, 0],
[W * math.cos(theta), W * math.sin(theta)],
[W * math.cos(theta) - H * math.sin(theta), W * math.sin(theta) + H * math.cos(theta)],
[-H * math.sin(theta), H * math.cos(theta)],
[0, 0]
])
# 计算水面
if theta <= enhanced_rect.alpha:
# 底部接触
water_points = np.array([
[0, 0],
[W * math.cos(theta), W * math.sin(theta)],
[W * math.cos(theta) - H * math.sin(theta), W * math.sin(theta) + H * math.cos(theta)],
[-H * math.sin(theta), H * math.cos(theta)]
])
elif theta < math.pi/2:
# 顶部接触
water_height = math.sqrt(2 * area * math.tan(theta))
water_points = np.array([
[0, 0],
[water_height * math.cos(theta), water_height * math.sin(theta)],
[-H * math.sin(theta), H * math.cos(theta)]
])
elif theta < math.pi:
# 倒置状态
water_points = np.array([
[0, 0],
[-H * math.sin(theta), H * math.cos(theta)],
[W * math.cos(theta) - H * math.sin(theta), W * math.sin(theta) + H * math.cos(theta)]
])
# 更新图形
container_line.set_data(corners[:, 0], corners[:, 1])
water_poly.set_xy(water_points)
# 更新文本
theta_text.set_text(f'Angle: {math.degrees(theta):.1f}°')
flow_text.set_text(f'Flow: {sim_results["flow_rate_actual"][frame]:.2f} m²/s (Target: {Q})')
return container_line, water_poly, pivot_point, theta_text, flow_text
# 创建动画
ani = FuncAnimation(fig, update, frames=len(sim_results['time']),
init_func=init, blit=True, interval=50)
plt.tight_layout()
plt.show()
在此代码基础下进行完善
最新发布