Doubt in Coupled Ode

2 次查看（过去 30 天）
Susmita Panda
Susmita Panda 2023-11-15
评论： Torsten 2023-11-16
Ran in:
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 个评论
Torsten
Torsten 2023-11-16
Either the solution you took from the book is wrong or your differential equations are wrong.

请先登录，再进行评论。

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

采纳的回答

Torsten
Torsten 2023-11-15
编辑：Torsten 2023-11-15
Ran in:
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
ans = 
sol.y
ans = 
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)
ans = 4×1
1.0e+03 * -2.1039 2.1039 0.1942 -0.1942
  1 个评论
Susmita Panda
Susmita Panda 2023-11-16
You are correct sir. I must check the values of constants.

请先登录，再进行评论。

更多回答（1 个）

Sam Chak
Sam Chak 2023-11-16
Ran in:
Hi @Susmita Panda
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)
ans = 4×1
1.0e+03 * 2.1039 -2.1039 -0.1942 0.1942
  2 个评论
Sam Chak
Sam Chak 2023-11-16
Ran in:
Hi @Susmita Panda
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
Susmita Panda
Susmita Panda 2023-11-16
编辑：Susmita Panda 2023-11-16
Thank you for the idea sir; however the problem in my case I am getting solution using analytical solution and not getting using ode45. I am also supposing that my constants are errorneous. I should provide the actual code and actual values of constant instead of demo code. Sir do check my revised question.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by