Doubt in Coupled Ode
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
I have written a code for coupled ode using ode45; however my results seems errornous when compared with analytical results:
%% Analytical
L=24;
A=9;
Iz=2.43;
J=21.18;
E=33.2e9;
EI=E*Iz;
rho=2.4e4;
volume=A*L;
m=(rho*A);
Fv=29.9e4;
theta=30*(pi/180);
R=L/theta;
v=40;
mu=0.2;
G=E/(2*(1+mu));
GJ=G*J;
g=10;
a1=(1/m)*(pi/L)^2*(EI*((pi/L)^2)+(GJ/(R^2)));
a2=(1/(m*R))*((pi/L)^2)*(EI+GJ);
b1=-(1/(rho*J))*((EI/R^2)+(GJ*((pi/L)^2)));
b2=-(1/(rho*J))*(1/R)*((pi/L)^2)*(EI+GJ);
wv=32.10; %% book formula provided in snip shot with constants
% factor_1=sqrt(((a1-b1)^2)+(4*a2*b2));
% wv=sqrt((a1+b1+factor_1)/2) ;% however I am getting different with my constants
g=10;
t=0.01:0.01:(L/v);
Sv=(pi*v)/(L*wv);
beta=(b1-(pi*v/L)^2)/(b1+a1-wv^2-(pi*v/L)^2);
xi=sin(pi*v*t/L)-(Sv*sin(wv*t));
u_deflect=-((2*Fv*g)/(m*L))*(1/wv^2)*(1/(1-Sv^2))*beta*xi;
figure();plot(t,u_deflect,'o-');title('Validation');
My code is:
%% Analysis using ode45
tspan_1=[0:0.001:0.6];%time range
y0_1=[0;0;0;0];%initial conditions f
[t1,y1]=ode45(@diffeqn11,tspan_1,y0_1);
%plot displacement
figure(1)
plot(t1, y1(:,1), 'r', 'LineWidth',2);title('Displacement of the beam');
hold on
plot(t1, y1(:,3), 'b', 'LineWidth',2);
%plot velocity_forced+free
figure(2)
plot(t1, y1(:,2), 'g', 'LineWidth',2);title('Velocity of the beam');
figure(3)
plot(t1, y1(:,4), 'r', 'LineWidth',2);
function f= diffeqn11(t,y)
L=24;
A=9;
Iz=2.43;
J=21.18;
E=33.2e9;
EI=E*Iz;
rho=2.4e4;
volume=A*L;
m=(rho*A);
Fv=29.9e4;
theta=30*(pi/180);
R=L/theta;
v=40;
mu=0.2;
G=E/(2*(1+mu));
GJ=G*J;
g=10;
a1=(1/m)*(pi/L)^2*(EI*((pi/L)^2)+(GJ/(R^2)));
a2=(1/(m*R))*((pi/L)^2)*(EI+GJ);
b1=-(1/(rho*J))*((EI/R^2)+(GJ*((pi/L)^2)));
b2=-(1/(rho*J))*(1/R)*((pi/L)^2)*(EI+GJ);
f=zeros(4,1);
f(1)=y(2);
f(2)=((2*Fv)/(m*L))*sin((pi*v*t)/L)-(a1*y(1))-(a2*y(3));
f(3)=y(4);
f(4)=-(b1*y(3))-(b2*y(1));
end
1 commentaire
Torsten
le 16 Nov 2023
Either the solution you took from the book is wrong or your differential equations are wrong.
Réponse acceptée
Torsten
le 15 Nov 2023
Modifié(e) : Torsten
le 15 Nov 2023
Your code is correct, but it seems there is a singularity near 0.32....
syms t x(t) y(t) a1 b1 a2 b2 Constant
eqn1 = diff(x,t,2)+a1*x+a2*y==Constant*sin(0.5*t);
eqn2 = diff(y,t,2)+b1*y+b2*x==0;
Dx = diff(x,t);
Dy = diff(y,t);
sol = dsolve([eqn1,eqn2],[x(0)==0,y(0)==0,Dx(0)==0,Dy(0)==0]);
sol.x
sol.y
figure(1)
fplot(subs(sol.x,[a1 a2 b1 b2 Constant],[120.7217,1.3587e+06,-4.4643e+06,-1.2327e+05,0.1154]),[0 0.6])
figure(2)
fplot(subs(sol.y,[a1 a2 b1 b2 Constant],[120.7217,1.3587e+06,-4.4643e+06,-1.2327e+05,0.1154]),[0 0.6])
Note that the eigenvalues of the characteristic polynomial for the system are tremendous:
a1 = 120.7217;
a2 = 1.3587e+06;
b1 = -4.4643e+06;
b2 = -1.2327e+05;
M = [0 0 1 0;0 0 0 1;-a1 -a2 0 0;-b2 -b1 0 0];
eig(M)
Plus de réponses (1)
Sam Chak
le 16 Nov 2023
Your 4th-order coupled system has two eigenvalues with positive real parts. Therefore, in theory, the state responses will diverge when the system is forced.
% parameters
a1 = 120.7217;
a2 = 1.3587e+06;
b1 = -4.4643e+06;
b2 = -1.2327e+05;
% state matrix
A = [ 0 1 0 0;
-a1 0 -a2 0;
0 0 0 1;
-b2 0 -b1 0];
% check eigenvalues
eig(A)
2 commentaires
Sam Chak
le 16 Nov 2023
Please note that in your beam-coupled system, there is typically no damping considered. However, in the real physical world, damping components exist. As a structural engineer, you can either incorporate other beams with stabilizing properties, or measure the deflections of existing beams and then counteract destabilizing vibration modes.
%% Analysis using ode45
tspan_1 = [0:0.001:4*pi]; % time range
y0_1 = [0; 0; 0; 0]; % initial conditions f
options = odeset('RelTol', 1e-8, 'AbsTol', 1e-8);
[t1, y1] = ode45(@diffeqn11, tspan_1, y0_1, options);
% plot displacement
figure(1)
plot(t1/pi, y1(:,1:2:3)); grid on,
title('Displacements of the beams'), xlabel('t/\pi')
legend('y_{1}', 'y_{3}', 'location', 'SE')
% A pair of Beam couples
function f = diffeqn11(t, y)
% parameters
a1 = 120.7217;
a2 = 1.3587e+06;
b1 = -4.4643e+06;
b2 = -1.2327e+05;
Constant= 0.1154;
% measure the deflections and velocities, and feed them back into the
% system in the same channel, where the sinusoidal force is injected.
k1 = 10562734.7949266;
k2 = 4596.24524333348;
k3 = 382536034.238245;
k4 = 181690.721406603;
K = [k1 k2 k3 k4];
f = zeros(4, 1);
u = Constant*sin(0.5*t) - K*y;
f(1) = 0*y(1) + 1*y(2) + 0*y(3) + 0*y(4);
f(2) = - a1*y(1) + 0*y(2) - a2*y(3) + 0*y(4) + u;
f(3) = 0*y(1) + 0*y(2) + 0*y(3) + 1*y(4);
f(4) = - b2*y(1) + 0*y(2) - b1*y(3) + 0*y(4);
end
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!