{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 蓖麻油粘滞系数计算 (斯托克斯定律)\n",
    "\n",
    "本项目用于计算落球法测量蓖麻油粘滞系数实验的数据处理与不确定度评定。\n",
    "\n",
    "---\n",
    "### 实验原理与计算公式\n",
    "\n",
    "#### 1. 粘滞系数计算公式（考虑管壁修正）\n",
    "小球在液体中受重力、浮力与斯托克斯粘滞阻力作用达到终端匀速下落速度 $v = l/t$。考虑到圆筒玻璃管壁的有限边界效应，引入管壁修正因子 $(1 + 2.4 d/D)$：\n",
    "$$\\eta = \\frac{(\\rho - \\rho_0) g d^2 t}{18 l \\left(1 + 2.4 \\frac{d}{D}\\right)}$$\n",
    "\n",
    "* $\\rho$：小球密度（钢球通常取 $7800\\,\\mathrm{kg/m^3}$）\n",
    "* $\\rho_0$：液体（蓖麻油）密度（室温下约 $960\\,\\mathrm{kg/m^3}$）\n",
    "* $g$：重力加速度（取 $9.80\\,\\mathrm{m/s^2}$）\n",
    "* $d$：小球直径\n",
    "* $D$：玻璃管内径\n",
    "* $l$：小球下落距离\n",
    "* $t$：下落平均时间\n",
    "\n",
    "#### 2. 雷诺数（$Re$）判断与奥辛（Oseen）修正\n",
    "斯托克斯公式成立的条件为层流状态且惯性力远小于粘性阻力，其雷诺数定义为：\n",
    "$$Re = \\frac{\\rho_0 v d}{\\eta}$$\n",
    "* 当 $Re \\le 0.1$ 时，斯托克斯定律完全适用，无需修正。\n",
    "* 当 $0.1 < Re \\le 1.0$ 时，流场出现微小扰动，采用奥辛公式进行一级修正：\n",
    "$$\\eta_{\\mathrm{corrected}} = \\frac{\\eta}{1 + \\frac{3}{16} Re}$$\n",
    "\n",
    "#### 3. 不确定度传递公式\n",
    "假定 $\\rho, \\rho_0, g$ 的误差可忽略不计，由误差合成法则，粘滞系数的相对不确定度传递公式为：\n",
    "$$u_r(\\eta) = \\frac{u_\\eta}{\\eta} = \\sqrt{ \\left(2 \\frac{u_d}{d}\\right)^2 + \\left(\\frac{u_t}{t}\\right)^2 + \\left(\\frac{u_l}{l}\\right)^2 + \\left(\\frac{u_D}{D}\\right)^2 }$$\n",
    "\n",
    "绝对不确定度为：\n",
    "$$u_\\eta = \\eta \\cdot u_r(\\eta)$$\n",
    "\n",
    "其中各物理量 $X$ 的不确定度按 A 类与 B 类合成：\n",
    "$$u_X = \\sqrt{u_{A,X}^2 + u_{B,X}^2}, \\quad u_{A,X} = \\frac{s_X}{\\sqrt{n}}, \\quad u_{B,X} = \\frac{\\Delta_{\\mathrm{inst}}}{\\sqrt{3}}$$"
   ],
   "id": "e5127401d3a8f15d"
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-09-22T07:54:14.369907400Z",
     "start_time": "2026-09-22T07:54:14.243997100Z"
    }
   },
   "source": [
    "import math\n",
    "from decimal import Decimal\n",
    "from python.utils import scientific_round, calculate_stats\n",
    "\n",
    "print(\"工具函数与模块加载完成。\")"
   ],
   "id": "ca5f57b71e825238",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "工具函数与模块加载完成。\n"
     ]
    }
   ],
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. 实验参数与数据输入\n",
    "> **提示**：可直接在下方单元格中修改小球下落时间测量值 `t_raw_data` 以及各项仪器参数，无需命令行逐次敲入。"
   ],
   "id": "4f13dab8eefdafa4"
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-09-22T07:54:26.075852500Z",
     "start_time": "2026-09-22T07:54:26.004272300Z"
    }
   },
   "source": [
    "# --- 物理常数与材料密度 ---\n",
    "G = Decimal(\"9.80\")          # 重力加速度 (m/s^2)\n",
    "RHO_BALL = Decimal(\"7800\")   # 小球(钢球)密度 (kg/m^3)\n",
    "RHO_OIL = Decimal(\"960\")     # 蓖麻油密度 (kg/m^3)\n",
    "\n",
    "# --- 几何与仪器参数 (单位: mm / cm) ---\n",
    "D_mean = Decimal(\"20.0\")      # 玻璃管直径 D (mm)\n",
    "d_mean = Decimal(\"1.0\")       # 小球直径 d (mm)\n",
    "l_mean = Decimal(\"20.0\")      # 下落距离 l (cm)\n",
    "\n",
    "# --- 仪器允差 (仪器误差) Δ_inst ---\n",
    "delta_D = Decimal(\"0\")         # D 的仪器误差 (mm)，如未单独评估可设为 0\n",
    "delta_d = Decimal(\"0\")         # d 的仪器误差 (mm)，如未单独评估可设为 0\n",
    "delta_l = Decimal(\"0\")         # l 的仪器误差 (mm)，如未单独评估可设为 0\n",
    "delta_t = Decimal(\"0.01\")      # 秒表仪器误差 (s)\n",
    "\n",
    "# --- 5次下落时间测量值 (s) ---\n",
    "# 请在此填入实际实验数据：\n",
    "t_raw_data = [Decimal(\"12.45\"), Decimal(\"12.38\"), Decimal(\"12.50\"), Decimal(\"12.42\"), Decimal(\"12.48\")]\n",
    "\n",
    "print(f\"参数设置完成:\")\n",
    "print(f\"- 玻璃管直径 D: {D_mean} mm\")\n",
    "print(f\"- 小球直径 d:   {d_mean} mm\")\n",
    "print(f\"- 下落距离 l:   {l_mean} cm\")\n",
    "print(f\"- 时间测量值 t:  {[float(x) for x in t_raw_data]} s\")"
   ],
   "id": "b222d81b56a956be",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "参数设置完成:\n",
      "- 玻璃管直径 D: 20.0 mm\n",
      "- 小球直径 d:   1.0 mm\n",
      "- 下落距离 l:   20.0 cm\n",
      "- 时间测量值 t:  [12.45, 12.38, 12.5, 12.42, 12.48] s\n"
     ]
    }
   ],
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. 测量数据统计与不确定度计算"
   ],
   "id": "b96c99013eae004b"
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-09-22T07:54:28.891288400Z",
     "start_time": "2026-09-22T07:54:28.810903200Z"
    }
   },
   "source": [
    "# 长度单位换算为 mm 以便统一不确定度计算\n",
    "l_mean_mm = l_mean * Decimal(\"10\")\n",
    "\n",
    "# 单次测量参数的 B 类不确定度 (u = Δ_inst / √3)\n",
    "D_u = delta_D / Decimal(str(math.sqrt(3)))\n",
    "d_u = delta_d / Decimal(str(math.sqrt(3)))\n",
    "l_u = delta_l / Decimal(str(math.sqrt(3)))\n",
    "\n",
    "# 5次时间测量的统计量 (平均值、合成不确定度、A类、B类)\n",
    "t_mean, t_u, t_ua, t_ub, t_s = calculate_stats(t_raw_data, delta_t)\n",
    "\n",
    "print(\"--- 时间统计分析 ---\")\n",
    "print(f\"时间平均值 t_mean : {t_mean:.4f} s\")\n",
    "print(f\"样本标准差 s(t)   : {t_s:.4f} s\")\n",
    "print(f\"A 类不确定度 u_A(t): {t_ua:.4f} s\")\n",
    "print(f\"B 类不确定度 u_B(t): {t_ub:.4f} s\")\n",
    "print(f\"合成不确定度 u(t)  : {t_u:.4f} s\")"
   ],
   "id": "b609644356d56555",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- 时间统计分析 ---\n",
      "时间平均值 t_mean : 12.4460 s\n",
      "样本标准差 s(t)   : 0.0477 s\n",
      "A 类不确定度 u_A(t): 0.0214 s\n",
      "B 类不确定度 u_B(t): 0.0058 s\n",
      "合成不确定度 u(t)  : 0.0221 s\n"
     ]
    }
   ],
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. 粘滞系数计算、雷诺数校验与修约输出"
   ],
   "id": "652642328be1ab03"
  },
  {
   "cell_type": "code",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2026-09-22T07:54:30.748035800Z",
     "start_time": "2026-09-22T07:54:30.683789300Z"
    }
   },
   "source": [
    "# 统一转换为国际单位制 (m, kg, s)\n",
    "D_m = D_mean / Decimal(\"1000\")\n",
    "d_m = d_mean / Decimal(\"1000\")\n",
    "l_m = l_mean_mm / Decimal(\"1000\")\n",
    "\n",
    "# --- 1. 斯托克斯粘滞系数计算 (带管壁修正) ---\n",
    "correction_factor = Decimal(\"1\") + Decimal(\"2.4\") * (d_m / D_m)\n",
    "numerator = (RHO_BALL - RHO_OIL) * G * (d_m**2) * t_mean\n",
    "denominator = Decimal(\"18\") * l_m * correction_factor\n",
    "eta_val = numerator / denominator\n",
    "\n",
    "# --- 2. 下落终端速度与雷诺数 Re 校验 ---\n",
    "v_val = l_m / t_mean\n",
    "Re = RHO_OIL * v_val * d_m / eta_val\n",
    "\n",
    "# --- 3. 雷诺数判断与 Oseen 修正 ---\n",
    "if Re > Decimal(\"0.1\"):\n",
    "    eta_corrected = eta_val / (Decimal(\"1\") + Decimal(\"3\") / Decimal(\"16\") * Re)\n",
    "    oseen_applied = True\n",
    "else:\n",
    "    eta_corrected = eta_val\n",
    "    oseen_applied = False\n",
    "\n",
    "# --- 4. 相对不确定度与绝对不确定度 ---\n",
    "# 各相对不确定度平方和项: (2*ud/d)^2 + (ut/t)^2 + (ul/l)^2 + (uD/D)^2\n",
    "term_d = (Decimal(\"2\") * d_u / d_mean) if d_mean != 0 else Decimal(\"0\")\n",
    "term_t = (t_u / t_mean) if t_mean != 0 else Decimal(\"0\")\n",
    "term_l = (l_u / l_mean_mm) if l_mean_mm != 0 else Decimal(\"0\")\n",
    "term_D = (D_u / D_mean) if D_mean != 0 else Decimal(\"0\")\n",
    "\n",
    "rel_u_sq = term_d**2 + term_t**2 + term_l**2 + term_D**2\n",
    "rel_u = Decimal(str(math.sqrt(float(rel_u_sq))))\n",
    "eta_u = eta_corrected * rel_u\n",
    "\n",
    "# --- 5. 按照科学修约规范对齐保留有效数字 ---\n",
    "eta_final, u_final = scientific_round(eta_corrected, eta_u)\n",
    "\n",
    "# --- 6. 实验报告格式化输出 ---\n",
    "print(\"=\" * 45)\n",
    "print(\"               实 验 计 算 结 果              \")\n",
    "print(\"=\" * 45)\n",
    "print(f\"下落平均时间 t    : {t_mean:.3f} ± {t_u:.3f} s\")\n",
    "print(f\"小球下落速度 v    : {v_val:.5f} m/s\")\n",
    "print(\"-\" * 45)\n",
    "print(f\"雷诺数 Re         : {Re:.4f}\")\n",
    "if not oseen_applied:\n",
    "    print(\"判定结论          : Re <= 0.1，满足斯托克斯层流条件，无需 Oseen 修正。\")\n",
    "else:\n",
    "    print(\"判定结论          : Re > 0.1，已自动应用 Oseen 修正公式。\")\n",
    "    print(f\"修正前粘滞系数 η0 : {eta_val:.4f} Pa·s\")\n",
    "print(\"-\" * 45)\n",
    "print(f\"相对不确定度 u_r  : {rel_u * 100:.2f}%\")\n",
    "print(f\"绝对不确定度 u_η  : {eta_u:.4f} Pa·s\")\n",
    "print(\"-\" * 45)\n",
    "print(f\"最终结果表示      : η = {eta_final} ± {u_final} Pa·s\")\n",
    "print(\"=\" * 45)"
   ],
   "id": "5abf7dd1e1769bf6",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=============================================\n",
      "               实 验 计 算 结 果              \n",
      "=============================================\n",
      "下落平均时间 t    : 12.446 ± 0.022 s\n",
      "小球下落速度 v    : 0.01607 m/s\n",
      "---------------------------------------------\n",
      "雷诺数 Re         : 0.0746\n",
      "判定结论          : Re <= 0.1，满足斯托克斯层流条件，无需 Oseen 修正。\n",
      "---------------------------------------------\n",
      "相对不确定度 u_r  : 0.18%\n",
      "绝对不确定度 u_η  : 0.0004 Pa·s\n",
      "---------------------------------------------\n",
      "最终结果表示      : η = 0.2069 ± 0.0004 Pa·s\n",
      "=============================================\n"
     ]
    }
   ],
   "execution_count": 4
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}