Force problem in FEA-SPH co-simulation

[Alex Lan]

Hi, I performed single wheel test with both rigid and flexible wheel referenced with demos. I got the wheel's force through FSI API. The functions are as follows:

std::vector<ChVector<>> ChSystemFsi::GetRigidBodiesForce() const {
    thrust::host_vector<Real3> forcesH = m_sysFSI->fsiGeneralData->rigid_FSI_ForcesD;
    std::vector<ChVector<>> forces;
    for (size_t i = 0; i < forcesH.size(); i++) {
forces.push_back(utils::ToChVector(forcesH[i]));
}
    return forces;
}

std::vector<ChVector<>> ChSystemFsi::GetRigidBodiesTorque() const {
thrust::host_vector<Real3> torquesH = m_sysFSI->fsiGeneralData ->rigid_FSI_TorquesD;
std::vector<ChVector<>> torques;
    for (size_t i = 0; i < torquesH.size(); i++) {
torques.push_back(utils::ToChVector(torquesH[i]));
}
    return torques;
}

std::vector<ChVector<>> ChSystemFsi::GetFlexBodiesForce() const {
thrust::host_vector<Real3> forcesH = m_sysFSI->fsiGeneralData->Flex_FSI_ForcesD;
std::vector<ChVector<>> forces;
    for (size_t i = 0; i < forcesH.size(); i++) {
forces.push_back(utils::ToChVector(forcesH[i]));
}
    return forces;
}

I got the force of the wheel like this:

void ReportForceTorqueFsi(ChSystemFsi& sysFSI, ChSystem& sysMBS, ChVector<>& m_force, ChVector<>& m_torque) {
    m_force = ChVector<>(0, 0, 0);
    m_torque = ChVector<>(0, 0, 0);

    if (flex_wheel) {
        auto nodes = sysFSI.GetFsiMesh()->GetNodes();
        std::vector<ChVector<>> node_force = sysFSI.GetFlexBodiesForce();

        auto wheel = sysMBS.Get_bodylist()[1];
        ChFrameMoving<> wheel_frame = *(wheel.get());

        for (int i = 0; i < nodes.size(); i++) {
            m_force += node_force[i];
            int index = nodes[i]->GetIndex();
            int m = index / (m_div_width + 1);
            int n = index % (m_div_width + 1);
            double phi = (CH_C_2PI * m) / m_div_circumference;
            double theta = -CH_C_PI * 2.0 / 3.0 + (CH_C_PI * 4.0 / 3.0 * n) / m_div_width;

            double x = (m_rim_radius + m_height * cos(theta)) * cos(phi);
            double y = m_height * sin(theta);
            double z = (m_rim_radius + m_height * cos(theta)) * sin(phi);
            ChVector<> pos = wheel_frame.TransformPointLocalToParent(ChVector<>(x, y, z)) - wheel_frame.TransformPointLocalToParent(ChVector<>(0, 0, 0));

            ChVector<> force = node_force[i];
            ChVector<> torque = Vcross(pos, force);
            m_torque += torque;
        }
    } else {
        auto force_list = sysFSI.GetRigidBodiesForce();
        auto torque_list = sysFSI.GetRigidBodiesTorque();
        for (int i = 0; i < force_list.size(); i++) {
m_force += force_list[i];
m_torque += torque_list[i];
        }
    }
}

Though the method to calculate torque may be not correct, the force of wheel is to add all the particle forces. For the two cases, flexible and rigid wheel, the wheels were both loaded 254N, like what the demos did. The reported force in vertical direction for rigid wheel is 254N correctly, but the force is always about 369N for flexible wheel.