Flapping Blade Problem

[Matt Tuman]

Hello Chrono Team,

I am attempting to model a rotor system with blades, but I am currently having issues capturing correct gyroscopic and inertial effects. To validate my approach, I am attempting to replicate a relatively simple system consisting of a “flapping blade”. The blade is rotating about its end and has one free DOF, \beta, as illustrated below. 

The only force in this system is the centrifugal force as illustrated. Integrating that force over the length of the blade results in the following equation of motion I_b \Omega^2 \beta + I_b \ddot{\beta} = 0, where I_b is the moment of inertia of the beam about its end. Thus, the natural frequency of the flapping mode is exactly equal to the rotation speed of the blade, \Omega. I am attempting to model this system in Project Chrono, but the modal solver continues to return a natural frequency of zero regardless of the rotation speed of the blade. I am assuming this issue can be fixed by properly initializing ChLoadBodyInertia, but I cannot find the solution. I have attached my code and corresponding output for reference. Please let me know if there is any additional information that would be helpful. Thanks!

Code:

int main(int argc, char* argv[]) {

    // Create a Chrono::Engine physical system

    ChSystemSMC sys;

    auto load_container = chrono_types::make_shared<ChLoadContainer>();

    sys.Add(load_container);

 

    // Define solver

    auto solver = chrono_types::make_shared<ChSolverPardisoMKL>();

    sys.SetSolver(solver);

 

    // Create ground

    auto ground = chrono_types::make_shared<ChBody>();

    ground->SetFixed(true);

    sys.Add(ground);

 

    // Create Rigid Body

    auto blade = chrono_types::make_shared<ChBodyEasyBox>(1., 0.1, 0.1, 0.0, true);

    blade->SetPos(ChVector3d(0.5, 0., 0.));

    blade->SetInertiaXX(ChVector3d(0., 0., 0.));

    blade->SetMass(0.);

    sys.Add(blade);

 

    // Add ChLoadBodyInertia to blade for gyroscopic effects

    auto gyro = chrono_types::make_shared<ChLoadBodyInertia>(blade, ChVector3d(0., 0., 0.), 1., ChVector3d(0.1, 1., 1.), ChVector3d(0., 0., 0.));

        // ChLoadBodyInertia(rigid body, CG offset, mass, InertiaXX, InertiaXY)

    load_container->Add(gyro);

 

    // Attach motor to spin blade about origin allowing it to "flap" about y axis

    auto motor = chrono_types::make_shared<ChLinkMotorRotationSpeed>();

    motor->Initialize(ground, blade, ChFrame<>(ChVector3d(0., 0., 0.)));

    std::shared_ptr<chrono::ChFunctionSetpoint> motor_speed = std::make_unique<chrono::ChFunctionSetpoint>();

    double speed = 10.; // [rad/s]

    motor_speed->SetSetpoint(speed, 0.);

    motor->SetSpeedFunction(motor_speed);  // Apply the angular velocity

    motor->SetSpindleConstraint(true, true, true, true, false); // allowing blade to rotate about y axis

    sys.AddLink(motor);

 

    // Perform static analysis

    ChStaticNonLinearRheonomicAnalysis static_solver;

    static_solver.SetIncrementalSteps(1);

    static_solver.SetMaxIterations(2);

    static_solver.SetResidualTolerance(1e-10);

    static_solver.SetVerbose(true);

    static_solver.SetAutomaticSpeedAndAccelerationComputation(true);

    sys.DoStaticAnalysis(static_solver);

 

    // Perform modal analysis

    chrono::ChMatrixDynamic<std::complex<double>> eigvects;

    chrono::ChVectorDynamic<std::complex<double>> eigvals;

    chrono::ChVectorDynamic<> freq;

    chrono::ChVectorDynamic<> damping_ratios;

    auto factorization = chrono_types::make_shared<ChSolverSparseComplexQR>();

    auto eig_solver = chrono_types::make_shared<ChUnsymGenEigenvalueSolverKrylovSchur>(factorization);

    ChModalSolverDamped modal_solver(4, 1e-5, true, false, eig_solver);

    modal_solver.Solve(sys.GetAssembly(), eigvects, eigvals, freq, damping_ratios);

 

    // Print out natural frequency

    std::cout << "Expected Natural Frequency: " << speed << "\n Predicted Natural Frequency: " << freq[0] * CH_2PI

              << std::endl;

 

    return 0;

}

Output:

[Chris VanDamme]

I am also interested in the answer to this question, is anyone from the Project Chrono development team able to assist with this?

Thanks,

Chris

--
You received this message because you are subscribed to the Google Groups "ProjectChrono" group.
To unsubscribe from this group and stop receiving emails from it, send an email to projectchron...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/projectchrono/c63d9283-7baa-4e13-8a1b-04d3d079bfc3n%40googlegroups.com.

[chao peng]

Hi guys,

The possible reason may be that in your code, the rotational speed of the rigid body 'blade' is zero although you have specified the speed via ChLinkMotorRotationSpeed.

Go to ChLoadBodyInertia::ComputeJacobian(ChState* state_x, ChStateDelta* state_w), and try to check whether m_jacobians->R and m_jacobians->K are nonzero. Since the COG offset in your test is given as zero, I would expect m_jacobians->K=0, but the gyroscopic damping matrix m_jacobians->R should be nonzero if the rotational speed is properly specified.

    // R gyroscopic damping matrix terms (6x6, split in 3x3 blocks for convenience)
    m_jacobians->R.setZero();
    if (this->use_inertial_damping_matrix_R) {
        //  Ri = [0, - m*[w~][c~] - m*[([w~]*c)~]  ; 0 , [w~][I] - [([I]*w)~]  ]
        m_jacobians->R.block(0, 3, 3, 3) = -this->mass * (wtilde * ctilde + ChStarMatrix33<>(wtilde * c_m));
        m_jacobians->R.block(3, 3, 3, 3) = wtilde * I - ChStarMatrix33<>(I * v_w);
    }

    // K inertial stiffness matrix terms (6x6, split in 3x3 blocks for convenience)
    m_jacobians->K.setZero();
    if (this->use_inertial_stiffness_matrix_K) {
        ChStarMatrix33<> atilde(a_w);  // [a~]
        // Ki_al = [0, -m*[([a~]c)~] -m*[([w~][w~]c)~] ; 0, m*[c~][xpp~] ]
        m_jacobians->K.block(0, 3, 3, 3) =
            -this->mass * ChStarMatrix33<>(atilde * c_m) - this->mass * ChStarMatrix33<>(wtilde * (wtilde * c_m));
        m_jacobians->K.block(3, 3, 3, 3) = this->mass * ctilde * ChStarMatrix33<>(a_x);
    }

static_solver.SetAutomaticSpeedAndAccelerationComputation(true) is designed to be able to derive the velocity and acceleration in the system wise. In this automatic derivation, the system stiffness matrix (K) is used to do a static analysis. But the inertial stiffness matrix (m_jacobians->K) is zero in your case, and the tangent stiffness matrix of constraints Kc is not activated (as the default setting), finally the system stiffness matrix should be zero, resulting in a failure in the automatic derivation for system velocity and acceleration.

I would suggest specifying the rotational speed of 'blade' explicitly by yourself through ChStaticNonLinearRheonomicAnalysis::IterationCallback(). See .\demo_FEA_beams_rotor.cpp for details in the code implementation.

By the way, in your test demo, I will also suggest activating Kc by motor->SetUseTangentStiffness(true). Kc might affect the eigenvalues.

Additionally, since the system size in your test demo is quite small, I am not sure whether the Krylov Schur eigensolver is able to find all eigenvalues correctly. To be safe, I would suggest sending the system matrices M K R Cq to MATLAB and solve the eigenvalues from MATLAB to do a comparison.

In the end, the attached paper investigated the influence of various stiffness matrices of a rotating rotor on the eigenvalues. Hope it is helpful for you.

Best regards,

Chao PENG

[Matt Tuman]

Hi Peng,

Thank you very much for providing all of this helpful feedback and information. The paper was also quite useful. Using the recommendations provided, I was able to get closer to the expected solution. However, to get the exact analytical solution, I had to set inertia values defined in ChBodyLoadInertia() to zero. That is, in my model, chrono can only solve for the expected solution when the blade is a point mass. Ideally, I would like to include inertial properties and still arrive at the expected solution. Do you have any ideas on how to resolve this issue? I plan on running this same case study with a FE blade to see if the same issue persists.

Thanks,

Matt

You received this message because you are subscribed to a topic in the Google Groups "ProjectChrono" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/projectchrono/MvCd3juUJnE/unsubscribe.
To unsubscribe from this group and all its topics, send an email to projectchron...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/projectchrono/d7fe3538-cdd6-4214-8c93-5d495afd00fan%40googlegroups.com.