水流冲重物问题

一个横截面积为$S$的水管向上冲水，冲出水管时的速度为$v_0$，往上放一个质量为$M$的重物，问这个重物能被水流冲多高？（假设一堆理想条件）（来自我弟问的一道中学物理竞赛题）

(1) 先分析稳态情况。指定时间$t$内，冲出水管时水的长度为$v_{0}t$，体积为$S{v}_{0}t$，质量为$m=\rho S{v}_{0}t$，重物距水管口高度为$H$，水流冲击到重物之前的速度为$v_1=\sqrt{v_0^2-2gH}$，水流动量为$m{v}_1$，重物冲量为$Mgt$，由$m{v}_1=Mgt$解得
$$\rho S{v}_0\sqrt{v_0^2-2gH}=Mg$$
即
$$H=\frac{1}{2g}\left(v_0^2-(\frac{Mg}{\rho S{v}_0}{)}^2\right)(1)$$
当$M>\frac{\rho S{v}_0^2}{g}$时$H=0$。
(2) 再分析动态情况。在微小时间$\text{d}t$内，设水流冲击前重物的速度为$v(t)$，冲击后重物的速度为$v(t)+\text{d}v(t)$，则
$$\begin{array}{l}M\text{d}v(t)=\rho S{v}_0\text{d}t\sqrt{v_0^2-2gh(t)}-Mg\text{d}t\\ \frac{\text{d}h(t)}{\text{d}t}=v(t)\end{array}$$
设$k_1=\rho S{v}_0$，得到微分方程如下
$$h^{\prime \prime }=\frac{\rho S{v}_0}{M}\sqrt{v_0^2-2gh}-g$$
如果重物太轻可能会出现重物冲出水流的情况，也就是说重物的高度$h$超出了水流的高度$\sqrt{v_0^2-2gh}$，此时的微分方程即为
$$h^{\prime \prime }=-g$$
$S=0.1$，$v_0=2$，$M=10$时重物的高度变化情况如下图

而这种情况下水流的高度只有0.2m。重物一直在运动停不下来，因此考虑一下加入空气阻力，
$$h^{\prime \prime }=-g-k{h}^{\prime }$$
其中$k$为空气阻力系数。重物的高度变化情况如下图，且改变$k$的值不影响稳定时重物的高度。将数据代入(1)式计算得到$H=0.1875$m。

不考虑空气阻力时的 simucpp 代码如下

#include "simucpp.hpp"
using namespace simucpp;
using namespace std;
int main()
{
    int err=0;
    constexpr double rho = 1000;  // 水的密度，1000kg/m^3
    constexpr double S = 0.1;  // 横截面积，m^2
    constexpr double v0 = 2;  // 水流初速度，m/s
    constexpr double g = 10;  // 重力加速度，m/s^2
    constexpr double M = 10;  // 重物质量，kg
    Simulator sim1(4);
    FMIntegrator(mv, &sim1);  // 输出重物速度
    FMIntegrator(mh, &sim1);  // 输出重物高度
    FMFcn(mfcn1, &sim1);
    FMOutput(mout1, &sim1);
    sim1.connect(mh, mfcn1);
    sim1.connect(mfcn1, mv);
    sim1.connect(mv, mh);
    sim1.connect(mh, mout1);
    mfcn1->Set_Function([](double u){
        double waterv = v0*v0-2*g*u;
        double ans = -g;
        return (waterv>0) ? ans+rho*S*v0*sqrt(waterv)/M : ans;
    });
    sim1.Set_SampleTime(0.01);
    sim1.Initialize();
    sim1.Simulate();
    return 0;
}

考虑空气阻力时的 simucpp 代码如下

#include "simucpp.hpp"
using namespace simucpp;
using namespace std;
int main()
{
    int err=0;
    constexpr double rho = 1000;  // 水的密度，1000kg/m^3
    constexpr double S = 0.1;  // 横截面积，m^2
    constexpr double v0 = 2;  // 水流初速度，m/s
    constexpr double g = 10;  // 重力加速度，m/s^2
    constexpr double M = 10;  // 重物质量，kg
    constexpr double kv = 10;  // 空气阻力系数
    Simulator sim1(4);
    FMIntegrator(mv, &sim1);  // 输出重物速度
    FMIntegrator(mh, &sim1);  // 输出重物高度
    FMFcnMISO(mfcn1, &sim1);
    FMOutput(mout1, &sim1);
    sim1.connect(mh, mfcn1);
    sim1.connect(mv, mfcn1);
    sim1.connect(mfcn1, mv);
    sim1.connect(mv, mh);
    sim1.connect(mh, mout1);
    mfcn1->Set_Function([](double *u){
        double waterv = v0*v0-2*g*u[0];
        double ans = -g-kv*u[1];
        return (waterv>0) ? ans+rho*S*v0*sqrt(waterv)/M : ans;
    });
    sim1.Set_SampleTime(0.01);
    sim1.Initialize();
    sim1.Simulate();
    PRINT_VALUE(mout1->Get_OutValue());
    return 0;
}