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
On Thu, Nov 21, 2024 at 6:35 PM chao peng <pengc1...@gmail.com> wrote:
> 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
>
> 在2024年11月22日星期五 UTC+8 01:27:07<cvand...@gmail.com> 写道：
>
>> 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
>>
>> On Thu, Nov 14, 2024 at 11:51 AM Matt Tuman <mtum...@gmail.com> wrote:
>>
>>> 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.
>>>
>>> [image: flapping_blade.png]
>>>
>>> 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[image: chrono_output.png]:
>>>
>>> --
>>> 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
>>> <https://groups.google.com/d/msgid/projectchrono/c63d9283-7baa-4e13-8a1b-04d3d079bfc3n%40googlegroups.com?utm_medium=email&utm_source=footer>
>>> .
>>>
>> --
> 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
> projectchrono+unsubscr...@googlegroups.com.
> To view this discussion visit
> https://groups.google.com/d/msgid/projectchrono/d7fe3538-cdd6-4214-8c93-5d495afd00fan%40googlegroups.com
> <https://groups.google.com/d/msgid/projectchrono/d7fe3538-cdd6-4214-8c93-5d495afd00fan%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>
--
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 projectchrono+unsubscr...@googlegroups.com.
To view this discussion visit
https://groups.google.com/d/msgid/projectchrono/CAO%2Bfy0QojfXbPNMnPDw3hQh0moo_Z7PdE%3D%3DrBccs0kNjZStrLw%40mail.gmail.com.
