SOR
145 lines
5.4 KiB
Python
145 lines
5.4 KiB
Python
import math
|
||
from decimal import Decimal
|
||
|
||
"""
|
||
Viscosity Coefficient Calculator (Stokes' Law)
|
||
蓖麻油粘滞系数计算器 (斯托克斯定律)
|
||
|
||
参考了杨氏模量代码的高精度与不确定度修约逻辑。
|
||
（注：由于实验提供的是同一距离20cm的5次重复时间测量，此处直接求平均计算，不适用逐差法。逐差法通常用于等间距的多组累积物理量测量。）
|
||
|
||
公式说明:
|
||
1. 粘滞系数计算公式 (考虑管壁修正):
|
||
η = (ρ - ρ0) * g * d^2 * t / [18 * l * (1 + 2.4 * d/D)]
|
||
其中:
|
||
ρ : 小球密度
|
||
ρ0 : 液体(蓖麻油)密度
|
||
d : 小球直径
|
||
D : 玻璃管内径
|
||
l : 下落距离
|
||
t : 下落平均时间
|
||
|
||
2. 雷诺数 (Re) 判断和修正公式:
|
||
Re = ρ0 * v * d / η
|
||
其中 v = l / t。
|
||
当 Re < 0.1 时，斯托克斯定律完全适用。
|
||
当 0.1 <= Re <= 1.0 时，可使用奥辛(Oseen)修正:
|
||
η_corrected = η / (1 + 3/16 * Re)
|
||
|
||
3. 不确定度传递公式 (相对不确定度):
|
||
假定 ρ, ρ0, g 的不确定度可忽略，则:
|
||
u_r(η) = u_η / η = sqrt( (2 * u_d / d)^2 + (u_t / t)^2 + (u_l / l)^2 + (u_D / D)^2 )
|
||
|
||
绝对不确定度:
|
||
u_η = η * u_r(η)
|
||
|
||
其中各个量 X 的不确定度 u_X = sqrt( u_A,X^2 + u_B,X^2 )
|
||
u_A,X = s_X / sqrt(n)
|
||
u_B,X = Δ_inst / sqrt(3)
|
||
"""
|
||
|
||
from utils import (
|
||
scientific_round,
|
||
calculate_stats,
|
||
input_data_group
|
||
)
|
||
|
||
def main():
|
||
print("========================================")
|
||
print(" 蓖麻油粘滞系数计算 (斯托克斯定律) ")
|
||
print("========================================")
|
||
|
||
# --- 1. 常量设置 ---
|
||
G = Decimal("9.80") # 重力加速度 m/s^2
|
||
RHO_BALL = Decimal("7800") # 钢球密度 kg/m^3 (默认为钢球，通常实验室用此值)
|
||
RHO_OIL = Decimal("960") # 蓖麻油密度 kg/m^3 (随温度变化，常温附近约 960)
|
||
|
||
print(f"重力加速度 g = {G} m/s^2")
|
||
print(f"小球密度 ρ = {RHO_BALL} kg/m^3")
|
||
print(f"蓖麻油密度 ρ0 = {RHO_OIL} kg/m^3")
|
||
|
||
# --- 2. 测量数据输入 ---
|
||
print("\n--- 仪器与实验参数设定 ---")
|
||
D_mean = Decimal("20.0")
|
||
d_mean = Decimal("1.0")
|
||
l_mean = Decimal("20.0") # 默认单位 cm，后续转化为 mm 或 m 计算
|
||
|
||
ans = input("使用默认参数 (玻璃管直径D=20mm, 小球直径d=1mm, 下落距离l=20cm) 吗？[Y/n]: ").strip().lower()
|
||
if ans == 'n':
|
||
D_mean = Decimal(input("玻璃管直径 D (mm): ").strip() or "20.0")
|
||
d_mean = Decimal(input("小球直径 d (mm): ").strip() or "1.0")
|
||
l_mean = Decimal(input("下落距离 l (cm): ").strip() or "20.0")
|
||
|
||
l_mean_mm = l_mean * Decimal("10") # 统一 cm 转为 mm
|
||
|
||
print("\n--- 输入仪器误差 Δ_inst ---")
|
||
delta_D = Decimal(input("D 的仪器误差 (mm, 默认 0): ").strip() or "0")
|
||
delta_d = Decimal(input("d 的仪器误差 (mm, 默认 0): ").strip() or "0")
|
||
delta_l = Decimal(input("l 的仪器误差 (mm, 默认 0): ").strip() or "0")
|
||
delta_t = Decimal(input("t(时间)的仪器误差 (s, 默认 0.01): ").strip() or "0.01")
|
||
|
||
# D, d, l 我们假设只进行了单次测量或作为已知常数处理，所以只算B类不确定度
|
||
D_u = delta_D / Decimal(str(math.sqrt(3)))
|
||
d_u = delta_d / Decimal(str(math.sqrt(3)))
|
||
l_u = delta_l / Decimal(str(math.sqrt(3)))
|
||
|
||
# 输入5个时间数据
|
||
t_vals = input_data_group("小球下落时间 t", 5, "s")
|
||
t_mean, t_u = calculate_stats(t_vals, delta_t)[:2]
|
||
|
||
# --- 3. 粘滞系数计算 ---
|
||
# 统一转化为国际标准单位 (m, kg, s)
|
||
D_m = D_mean / Decimal("1000")
|
||
d_m = d_mean / Decimal("1000")
|
||
l_m = l_mean_mm / Decimal("1000")
|
||
|
||
# η = (ρ - ρ0) * g * d^2 * t / [18 * l * (1 + 2.4 * d/D)]
|
||
correction_factor = Decimal("1") + Decimal("2.4") * (d_m / D_m)
|
||
numerator = (RHO_BALL - RHO_OIL) * G * (d_m**2) * t_mean
|
||
denominator = Decimal("18") * l_m * correction_factor
|
||
eta_val = numerator / denominator
|
||
|
||
# --- 4. 雷诺数判断与修正 ---
|
||
# v = l / t
|
||
v_val = l_m / t_mean
|
||
# Re = ρ0 * v * d / η
|
||
Re = RHO_OIL * v_val * d_m / eta_val
|
||
|
||
eta_corrected = eta_val
|
||
if Re > Decimal("0.1"):
|
||
# Oseen 修正
|
||
eta_corrected = eta_val / (Decimal("1") + Decimal("3") / Decimal("16") * Re)
|
||
|
||
# --- 5. 不确定度传播 ---
|
||
# 相对不确定度平方和
|
||
rel_u_sq = (Decimal("2") * d_u / d_mean)**2 + (t_u / t_mean)**2 + (l_u / l_mean_mm)**2 + (D_u / D_mean)**2
|
||
rel_u = Decimal(str(math.sqrt(float(rel_u_sq))))
|
||
|
||
# 绝对不确定度
|
||
eta_u = eta_corrected * rel_u
|
||
|
||
# --- 6. 修约与输出 ---
|
||
eta_final, u_final = scientific_round(eta_corrected, eta_u)
|
||
|
||
print("\n" + "="*40)
|
||
print(" 计算结果")
|
||
print("="*40)
|
||
print(f"下落时间 t = {t_mean:.3f} ± {t_u:.3f} s")
|
||
print(f"下落速度 v = {v_val:.5f} m/s")
|
||
print("-" * 40)
|
||
print(f"雷诺数 Re = {Re:.4f}")
|
||
if Re <= Decimal("0.1"):
|
||
print("Re <= 0.1，满足斯托克斯定律条件，无需 Oseen 修正。")
|
||
else:
|
||
print("Re > 0.1，已应用 Oseen 修正。")
|
||
print(f"修正前 η = {eta_val:.4f} Pa·s")
|
||
print("-" * 40)
|
||
print(f"相对不确定度 u_r = {rel_u*100:.2f}%")
|
||
print(f"绝对不确定度 u_η = {eta_u:.4f} Pa·s")
|
||
print("-" * 40)
|
||
print(f"最终结果: η = {eta_final} ± {u_final} Pa·s")
|
||
print("="*40)
|
||
|
||
if __name__ == "__main__":
|
||
main()
|