指定物体运动或变形

SPHinXsys通过MotionConstraint类模板来实现对物体的（刚体）运动或（弹性体）变形控制。写法很简单，提供每个粒子的坐标、速度、加速度即可。坐标和速度是必需的，加速度非必需，但是有的动力学会用到加速度，最好一并提供。

MotionConstraint

/**
 * @class	MotionConstraint
 * @brief	Base class for constraining with prescribed motion.
 * 			Exact motion function will be defined in derive class.
 * 			Note that, we do not impose acceleration, so that this constraint
 * 			must be imposed after updating particle velocity by forces
 * 			and before updating particle position.
 */
template <class DynamicsIdentifier>
class MotionConstraint : public BaseLocalDynamics<DynamicsIdentifier>
{
  public:
    explicit MotionConstraint(DynamicsIdentifier &identifier)
        : BaseLocalDynamics<DynamicsIdentifier>(identifier),
          pos_(this->particles_->template getVariableDataByName<Vecd>("Position")),
          pos0_(this->particles_->template registerStateVariableDataFrom<Vecd>("InitialPosition", "Position")),
          vel_(this->particles_->template registerStateVariableData<Vecd>("Velocity")) {};

    virtual ~MotionConstraint() {};

  protected:
    Vecd *pos_, *pos0_, *vel_;
};

这是个比较简单的类模板。DynamicsIdentifier一般会设定为SPHBody或者BodyPartByParticle这些有LoopRange成员函数的类。BaseLocalDynamics<DynamicsIdentifier>是基类，提供粒子属性的访问接口。在构造时，会从body拿到所属粒子的当前位置pos_、初始位置pos0_和速度vel_。

using BodyPartMotionConstraint = MotionConstraint<BodyPartByParticle>;

如果DynamicsIdentifier是BodyPartByParticle，可以使用别名MotionConstraint<BodyPartByParticle>。

FixConstraint和FixedInAxisDirection是两个派生类。前者在所有方向上固定住指定的body或body part；后者在给定方向（x/y/z）上固定住指定的body或body part。

案例

这里提供两个案例，一个是2d_examples/test_2d_owsc案例对造波器的平动控制，另一个是我自己写的蠕动泵壁面变形控制。

造波器的平动

class WaveMaking : public BodyPartMotionConstraint
{
    Real model_scale_;
    Real gravity_;
    Real water_depth_;
    Real wave_height_;
    Real wave_period_;
    Real wave_freq_;
    Real wave_stroke_;

    Vecd getDisplacement(const Real &time)
    {
        Vecd displacement{Vecd::Zero()};
        displacement[0] = 0.5 * wave_stroke_ * sin(wave_freq_ * time);
        return displacement;
    }

    Vec2d getVelocity(const Real &time)
    {
        Vec2d velocity{Vecd::Zero()};
        velocity[0] = 0.5 * wave_stroke_ * wave_freq_ * cos(wave_freq_ * time);
        return velocity;
    }

    Vec2d getAcceleration(const Real &time)
    {
        Vec2d acceleration{Vecd::Zero()};
        acceleration[0] = -0.5 * wave_stroke_ * wave_freq_ * wave_freq_ * sin(wave_freq_ * time);
        return acceleration;
    }

    void computeWaveStrokeAndFrequency()
    {
		// 计算wave_stroke_与wave_freq_
    }

  public:
    WaveMaking(BodyPartByParticle &body_part)
        : BodyPartMotionConstraint(body_part),
          model_scale_(25.0), gravity_(gravity_g), water_depth_(Water_H), wave_height_(5.0),
          wave_period_(10.0),
          acc_(particles_->registerStateVariableData<Vecd>("Acceleration")),
          physical_time_(sph_system_.getSystemVariableDataByName<Real>("PhysicalTime"))
    {
        computeWaveStrokeAndFrequency();
    }

    void update(size_t index_i, Real dt = 0.0)
    {
        Real time = *physical_time_;
        pos_[index_i] = pos0_[index_i] + getDisplacement(time);
        vel_[index_i] = getVelocity(time);
        acc_[index_i] = getAcceleration(time);
    };

  protected:
    Vecd *acc_;
    Real *physical_time_;
};

这里最重要的就是构造函数和update函数，其余都是辅助函数。因为造波器是wall_boundary的一部分，所以让WaveMaking继承于BodyPartMotionConstraint。额外指定了acc_和physical_time_两个属性，用于获取当前时间和指定加速度。update中，首先获取了当前的物理时间，然后为当前粒子赋予计算得到的坐标、速度和加速度。从辅助函数中，不难看出位移就是一个简单的正弦函数，速度和加速度分别是一阶和二阶导数。

按照以下方式使用。因为WaveMaking只需要update，所以使用SimpleDynamics即可。

int main()
{
    ...
    BodyRegionByParticle wave_maker(wall_boundary, makeShared<MultiPolygonShape>(createWaveMakerShape()));
    SimpleDynamics<WaveMaking> wave_making(wave_maker);
    ...
    while(...){
        ...
        wave_making.exec(dt);
        ...
    }
    ...
}

蠕动泵的变形

主要参考了COMSOL的蠕动泵案例（只有单层壁面，没有流体）：

公式差不多是照搬过来的，只不过我改了改参数：

参数：

//	Basic geometry parameters and numerical setup.
//----------------------------------------------------------------------
const Real l_si = 0.04;
const Real r_si = 0.002;
const Real resolution_ref = 2e-4;
const Real wall_thickness = 3e-4;
const int n_ring = int(l_si/resolution_ref);
const Real l_buf = r_si * 0.5;
// peristalsis params
const Real t_on = 1;
const Real t_off = 3;
const Real tc = 7;
const Real width2 = 0.004*0.004;
const Real x0 = 0.001;
const Real v0 = 0.008;
const Real strain_max = 0.4;
const Real cycle = 4;

壁面几何采用的是2d_examples/test_2d_channel_fluid_shell的壁面几何（参数改了）：

namespace SPH
{
class WallBoundary;
template <>
class ParticleGenerator<SurfaceParticles, WallBoundary> : public ParticleGenerator<SurfaceParticles>
{
  public:
    explicit ParticleGenerator(SPHBody &sph_body, SurfaceParticles &surface_particles)
        : ParticleGenerator<SurfaceParticles>(sph_body, surface_particles) {}
    void prepareGeometricData() override
    {
        for (int i = 0; i < n_ring; i++)
        {
            Real x = (Real(i) + 0.5) * resolution_ref;
            // upper wall
            Real y1 = r_si + 0.5 * resolution_ref;
            addPositionAndVolumetricMeasure(Vecd(x, y1), resolution_ref);
            Vec2d normal_direction_1 = Vec2d(0, 1.0);
            addSurfaceProperties(normal_direction_1, wall_thickness);
            // lower wall
            Real y2 = -r_si - 0.5 * resolution_ref;
            addPositionAndVolumetricMeasure(Vecd(x, y2), resolution_ref);
            Vec2d normal_direction_2 = Vec2d(0, -1.0);
            addSurfaceProperties(normal_direction_2, wall_thickness);
        }
    }
};
} // namespace SPH

变形控制类。因为蠕动泵这个变形函数太复杂了，难求导，所以我没有求解析的导数，而是用中心差分计算了近似导数。另外，我想控制整个壁面的变形，因此无需用body part，而是直接用SPHBody。所以这个类继承于MotionConstraint<SPHBody>。

class PeristalticContraction : public MotionConstraint<SPHBody>
{
  private:
    Vecd *acc_;
    Real *physical_time_;
    Real delta_t_; // for numerical differentiation

    Real getDisplacement(const Real &time, const Real &x)
    {
        Real t_local = std::fmod(time, cycle);
        Real dist_peak = x - (x0 + v0 * t_local);
        Real tmp = flc2hs(t_off - t_local, 1) * flc2hs(t_local - t_on, 1) * 
                   std::exp(-dist_peak*dist_peak / (2 * width2));
        return strain_max * r_si * tmp;
    }

    Real getVelocity(const Real &time, const Real &x)
    {
        Real disp_prev = getDisplacement(time - delta_t_, x);
        Real disp_next = getDisplacement(time + delta_t_, x);
        return (disp_next - disp_prev) / (2 * delta_t_);
    }

    Real getAcceleration(const Real &time, const Real &x)
    {
        Real vel_prev = getVelocity(time - delta_t_, x);
        Real vel_next = getVelocity(time + delta_t_, x);
        return (vel_next - vel_prev) / (2 * delta_t_);
    }

    static Real flc2hs(Real x, Real d)
    {
        if (x <= -d) return 0;
        if (x >= d) return 1;
        return  0.5 + 0.9375 * x - 0.625 * x*x*x + 0.1875 * x*x*x*x*x;
    }

  public:
    PeristalticContraction(SPHBody &body, Real delta_t)
        : MotionConstraint<SPHBody>(body),
          acc_(particles_->registerStateVariableData<Vecd>("Acceleration")),
          physical_time_(sph_system_.getSystemVariableDataByName<Real>("PhysicalTime")),
          delta_t_(delta_t) {}

    void update(size_t index_i, Real dt = 0.0)
    {
        Real time = *physical_time_;
        Real direction = pos0_[index_i][1] < 0 ? 1 : -1;
        pos_[index_i][1] = pos0_[index_i][1] + getDisplacement(time, pos_[index_i][0]) * direction;
        vel_[index_i][1] = getVelocity(time, pos_[index_i][0]) * direction;
        acc_[index_i][1] = getAcceleration(time, pos_[index_i][0]) * direction;
    }
};

当前案例只有固体变形，没有流体。如果有流体，并且没有实现薄壳动力学的话，还要在这里考虑加入对法线方向的指定；如果实现了薄壳动力学，可以用UpdateShellNormalDirection更新法线方向。

变形类的使用：

int main()
{
    /* Domain bounds of the system. */
    BoundingBox system_domain_bounds(Vec2d(0, -r_si * 1.5 - 0.5 * resolution_ref),
                                     Vec2d(l_si, r_si * 1.5 + 0.5 * resolution_ref));
    SPHSystem sph_system(system_domain_bounds, resolution_ref);
    //	Creating body, materials and particles.
    SolidBody wall_boundary(sph_system, makeShared<BinaryShapes>("Wall"));
    wall_boundary.defineMaterial<Solid>();
    wall_boundary.generateParticles<SurfaceParticles, WallBoundary>();
    // movement 
    SimpleDynamics<PeristalticContraction> peristaltic_contract(wall_boundary, 1e-5);
    // IO
    BodyStatesRecordingToVtp body_states_recording(sph_system);
    body_states_recording.addToWrite<Vecd>(wall_boundary, "Velocity");
    body_states_recording.addToWrite<Vecd>(wall_boundary, "Acceleration");
    body_states_recording.writeToFile();
	// main loop
    Real &physical_time = *sph_system.getSystemVariableDataByName<Real>("PhysicalTime");
    Real end_time = 7;
    Real dt = 0.02;
    while (physical_time < end_time)
    {
        peristaltic_contract.exec();
        physical_time += dt;
        std::cout << std::fixed << std::setprecision(9) << physical_time << std::endl;
        body_states_recording.writeToFile();
    }

    return 0;
}

下图展示了两个时刻的壁面Y方向速度：