Divergent solution using ode45

9 次查看（过去 30 天）
Abhik Saha
Abhik Saha 2022-9-9
评论： Sam Chak 2024-7-12
I have a code (see below) that solve the second order differential equation. As a output it gives the position and velocity as a function of time. For some values of omega 0.48,0.5,0.85 it gives the diverging solution of position. I guess this is due to some error but I can not find the error. How to overcome this errors for above omegas. Please help regarding this. I am new to MATLAB.
%% MAIN PROGRAM %%
t=linspace(0,1000,5000);
y0=[1 0];
omega=0.85;
[tsol, ysol]=ode45(@(t,y0) firstodefun4(t,y0,omega), t, y0, omega);
position=ysol(:,1);
velocity=ysol(:,2);
%% FUNCTION DEFINITION%%
function dy=firstodefun4(t,y0,omega)
F=1;gamma=0.01;omega0=1;
kappa=0.15;
dy=zeros(2,1);
dy(1)=y0(2);
dy(2)=2*F*sin(omega*t)-2*gamma*y0(2)-omega0^2*(1+4*kappa*sin(2*omega*t))*y0(1);
end

请先登录，再回答此问题。

回答（2 个）

Torsten
Torsten 2022-9-9
编辑：Torsten 2022-9-9
Ran in:
We don't know the reason - technically, your code is correct. Although I would change it as shown below.
What makes you think that the solution you get is incorrect ? Do you have any indication that it should be bounded as t -> 00 ? For me, the behaviour looks very regular.
%% MAIN PROGRAM %%
t=linspace(0,50,5000);
y0=[1 0];
omega=0.85;
F=1;
gamma=0.01;
omega0=1;
kappa=0.15;
[tsol, ysol]=ode45(@(t,y0) firstodefun4(t,y0,omega,F,gamma,omega0,kappa), t, y0);
position=ysol(:,1);
velocity=ysol(:,2);
plot(tsol,position)
%% FUNCTION DEFINITION%%
function dy=firstodefun4(t,y0,omega,F,gamma,omega0,kappa)
dy=zeros(2,1);
dy(1)=y0(2);
dy(2)=2*F*sin(omega*t)-2*gamma*y0(2)-omega0^2*(1+4*kappa*sin(2*omega*t))*y0(1);
end
  4 个评论
显示 2更早的评论
Abhik Saha
Abhik Saha 2022-9-9
I have checked with RelTol and Abstol but it does not change much I have also check with smaller time steps but again the result remains same. If possible can you try with RADAU5 from Hairer-Wanner because I do not know how to use this. HJust to check whether it gives divergent result or not. Apart from that I have not tried with kappa=0. But I have tried with other ode integrator it does not change much.
Jano
Jano 2022-10-14
Is it possible for this to accept a range of value for omega, for example omega = [0.85:5].

请先登录，再进行评论。

Sam Chak
Sam Chak 2022-9-12
Ran in:
Hi @Abhik Saha
If you reduce the value for , then the system should be stable.
For more info, please check out Floquet theory.
% MAIN PROGRAM %
n = 400;
t = linspace(0, n*100, n*10000+1);
y0 = [1 0];
omega = 0.85;
F = 1;
gamma = 0.01;
omega0 = 1;
kappa = 0.1469;
[tsol, ysol] = ode45(@(t, y) firstodefun4(t, y, omega, F, gamma, omega0, kappa), t, y0);
position = ysol(:, 1);
velocity = ysol(:, 2);
plot(tsol, position)
% ODE
function dy = firstodefun4(t, y, omega, F, gamma, omega0, kappa)
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = 2*F*sin(omega*t) - 2*gamma*y(2) - omega0^2*(1 + 4*kappa*sin(2*omega*t))*y(1);
end
  5 个评论
显示 3更早的评论
Nikoo
Nikoo 2024-7-12
编辑：Sam Chak 2024-7-12
Ran in:
@Sam Chak
I was implementing the unstability and divergence boundries of a system with 2dof in MATLAB, but my matrices have periodic coefficient.
clc
clear all ;
% Define parameters
N = 2 ; %Number of blades
I_thetaoverI_b = 2 ; % Moment of inertia pitch axis over I_b
I_psioverI_b = 2 ; % Moment of inertia yaw axis over I_b
C_thetaoverI_b = 0.00; % Damping coefficient over I_b
C_psioverI_b = 0.00; % Damping coefficient over I_b
h = 0.3; % rotor mast height, wing tip spar to rotor hub
hoverR = 0.34;
R = h / hoverR;
gamma = 4; % lock number
V = 325 ; % the rotor forward velocity [knots]
Omega = V/R; % the rotor rotational speed [RPM]
freq_1_over_Omega = 1 / Omega;
%the flap moment aerodynamic coefficients for large V
M_b = -(1/10)*V;
M_u = 1/6;
%the propeller aerodynamic coefficients
H_u = V/2;
%%%%%%%%%%%the flap moment aerodynamic coefficients for small V
%M_b = -1*(1 + V^2)/8 ;
%M_u = V/4;
%the propeller aerodynamic coefficients
%H_u = (V^2/2)*log(2/V);
f_pitch= 0.01:5:140;
f_yaw= 0.01:5:140;
phi_steps = linspace(0, pi, 50); % Evaluation points from 0 to pi
divergence_map = [];
Rdivergence_map = [];
unstable = [];
for i = 1:length(f_pitch)
for j = 1:length(f_yaw)
for phi = phi_steps
% Calculate stiffness for the current frequency
w_omega_pitch = 2*pi*f_pitch(i);
w_omega_yaw = 2*pi*f_yaw(j);
K_psi = (w_omega_pitch^2) * I_psioverI_b;
K_theta = (w_omega_yaw^2) * I_thetaoverI_b;
% Define inertia matrix [M]
M_matrix = [I_thetaoverI_b + 1 + cos(2*phi), -sin(2*phi);
-sin(2*phi), I_psioverI_b + 1 - cos(2*phi)];
% Define damping matrix [D]
D11 = h^2*gamma*H_u*(1 - cos(2*phi)) - gamma*M_b*(1 + cos(2*phi)) - (2 + 2*h*gamma*M_u)*sin(2*phi);
D12 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) - 2*(1 + cos(2*phi)) - 2*h*gamma*M_u*cos(2*phi);
D21 = h^2*gamma*H_u*sin(2*phi) + gamma*M_b*sin(2*phi) + 2*(1 - cos(2*phi)) - 2*h*gamma*M_u*cos(2*phi);
D22 = h^2*gamma*H_u*(1 + cos(2*phi)) - gamma*M_b*(1 - cos(2*phi)) + (2 + 2*h*gamma*M_u)*sin(2*phi);
D_matrix = [D11, D12;
D21, D22];
% Define stiffness matrix [K]
K11 = K_theta - h*gamma*V*H_u*(1 - cos(2*phi)) + gamma*V*M_u*sin(2*phi);
K12 = -h*V*gamma*H_u*sin(2*phi) + gamma*V*M_u*(1 + cos(2*phi));
K21 = -h*gamma*V*H_u*sin(2*phi) - gamma*V*M_u*(1 - cos(2*phi));
K22 = K_psi - h*gamma*V*H_u*(1 + cos(2*phi)) - gamma*V*M_u*sin(2*phi);
K_matrix = [K11, K12;
K21, K22];
A_top = [zeros(2, 2), eye(2)];
A_bottom = [-inv(M_matrix) * K_matrix, -inv(M_matrix) * D_matrix];
A = [A_top; A_bottom];
eigenvalues = eig(A);
% Stability condition
% Flutter
if any(real(eigenvalues) > 0)
unstable = [unstable; K_psi, K_theta];
end
% Divergence condition
if det(K_matrix) < 0
divergence_map = [divergence_map; K_psi, K_theta];
end
% 1/Ω *(Divergence) proximity check
for ev = eigenvalues'
if abs(ev - freq_1_over_Omega) < 1e-2
Rdivergence_map = [Rdivergence_map; K_psi, K_theta];
end
end
end
end
end
% Plot the Flutter and divergence maps
figure;
hold on;
scatter(unstable(:,1), unstable(:,2), 'filled');
scatter(divergence_map(:,1), divergence_map(:,2), 'filled', 'r');
scatter(Rdivergence_map(:, 1), Rdivergence_map(:, 2), 'filled', 'g');
xlabel('K\_psi');
ylabel('K\_theta');
title('Whirl Flutter Diagram');
legend('Flutter area', 'Divergence area', ' 1/Ω Divergence area');
hold off;
Could you take a look at my approach? The final plot isn't implemented correctly.
Sam Chak
Sam Chak 2024-7-12
Hi @Nikoo
Some engineering stuff here as I can see the eigenvalues of a two coupled 2nd-order system. I would suggest you to post as New Question to attract more attentions from the Aeroelastic flutter experts. If possible, attach a sketch of the expected result or plot.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by