Hi Nikko,
First, let me say that you did a very good job with the code. Now, going back to your comment about “final plot didn't match the reference”, what does your reference plot looks like, please elaborate little further by providing more details in context to that why your final plot does not match the reference plot? Upon reviewing your code, it appears that you're attempting to create a whirl flutter diagram based on these calculations by identifying regions of instability such as flutter and divergence areas by examining the eigenvalues of the system. Normally, the visualization of these regions is done through scatter plots. Without looking at your reference plot, I did analyze your code step by step and this is what I came across, the parameters you have defined seem fine and do not require any changes, the nested loops you have used to iterate over f_pitch, f_yaw, and phi_steps are correct. They allow you to calculate the stability for different frequencies and evaluation points, calculation of the inertia matrix seems correct and does not require any changes, calculation of the damping matrix also seems correct and does not require any changes, calculation of the stiffness matrix seems correct and does not require any changes, assembly of the state-space matrix seems correct and does not require any changes, calculation of the eigenvalues and the stability analysis seem fine.
However, there is one improvement you can make. Instead of using any(real(eigenvalues) > 0) to check for flutter, you can use the max(real(eigenvalues)) to find the maximum real part of the eigenvalues. If the maximum real part is positive, it indicates flutter.
The calculation of the divergence condition seems correct and does not require any changes, proximity check for 1/Ω *(Divergence) seems fine and does not require any changes and plotting code seems fine and does not require any changes.
Now, to address your query regarding use of Floquet technique and how to modify the current code, you have to follow these steps,
In your provided code, the frequency range is defined using the variables f_pitch and f_yaw, you need to define the time period of the system instead. Let's call it T. You can calculate T using the formula T = 1 / (2 * pi * f), where f is the frequency of the system. Modify the code to calculate T for both pitch and yaw frequencies. Modify the loop for evaluating the system at different time points, the loop iterates over the variables f_pitch and f_yaw to evaluate the system at different frequencies, you need to modify the loop to iterate over different time points within one time period T. Let's call the time variable t. Modify the loop to iterate over t from 0 to T with a certain step size, and within the modified loop, calculate the stability matrix A at each time point t. The stability matrix is defined as the derivative of the system's state vector with respect to time. Make sure that in your provided code, the stability matrix is calculated using the variables M_matrix, D_matrix, and K_matrix and update these matrices based on the modified code. After calculating the stability matrix A at each time point, calculate its eigenvalues using the eig() function in MATLAB. The eigenvalues represent the stability characteristics of the system. If any of the eigenvalues have a positive real part, it indicates instability. Based on the calculated eigenvalues, analyze the stability and divergence conditions of the system. In your code, the variables unstable and divergence_map are used to store the values of K_psi and K_theta for unstable and divergent conditions, respectively, modify these variables to store the values of t and K_psi or K_theta for the corresponding conditions. Finally, plot the Flutter and divergence maps using the scatter() function in MATLAB. Bear in mind that the x-axis represents K_psi, and the y-axis represents K_theta by using the modified variables unstable, divergence_map, and Rdivergence_map.
Please let me know if you have any further questions.
