Nondimentionalize the 2 DOF equations
11 views (last 30 days)
Show older comments
Hello all,
I hope you are doing well. I asked a question regarding dimensionalizing the motion equations in 2DOFs. Due to my mistakes, the description is not very clear. As suggested, I post a single question for the same question.
I would like to present the system I want to establish in MATLAB firstly. The absorber consists of a mass, three-spring system in the left fig, and the dashpot.
For the angular velocity , damping ratioξ, mass ratio μ, and the stiffness ratio α, they are easy to determine. The length of the oblique spring in my code is set to be 2.5m. But it can be various physically.
I am not sure the approach to achieve the nondimensionlized equation is corret or not. Also, I would like to know whether I could use the ratio value only. Specifically, I do not need to set the value.
My codes are also attached. Thank you for your help.
Best wishes,
Yu
clc;
clear all;
tspan = 0:0.25:200;
% Units:
% stiffness: N/m, mass:kg, damping(c_1):Ns/m
% displacement(X):m
% alpha (stiffness ratio):k_o/k_v
% gamma (dimension ratio):a/L (structure length)
% f_opt (frequency ratio): f_1/f_2
% k_v the stiffness of the vertical spring
% k_o the stiffness of the oblique spring
X0 = [0 0 0 0];
%parameters-------
global alpha gamma mu Omega_1 Omega_2 L_0
alpha = 1 % stiffness ratio k_o/k_v
gamma = 2*alpha/(2*alpha+1)
L_0 = 2.5
a = L_0*gamma;
h = sqrt(L_0^2-a^2);
mu = 0.02; % mass ratio
f_opt = 1/(1+mu); % frequency ratio
xi_2 = sqrt(3*mu/8*(1+mu));
Omega_1 = 0.188; % 0.03 Hz angular velocity of the primary structure
Omega_2 = Omega_1*f_opt; % angular velocity of the damper
xi_1 = 0.01; % damping ratio of primary structure
%matrix--------
M = [(1+mu)/(Omega_1^2*L_0) mu/(Omega_1^2*L_0);
1/(Omega_2^2*L_0) 1/(Omega_2^2*L_0)];
C = [2*xi_1/(Omega_1*L_0) 0;
0 2*xi_2/(Omega_2*L_0)];
K = [1/L_0 0;
0 1/L_0];
% state space model-------------------
% M: mass matrix, K:stiffness matrix, C:damping matrix
O = zeros(2,2);
I = eye(2);
A = [O I; -inv(M)*K -inv(M)*C];
B = [O; inv(M)];
E = [I O];
D = zeros(2,2);
% solve the equations-------------------
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t,X] = ode45(@(t,X) QZSdamper(t,A,B,X),tspan,X0,options);
x = X(:,1:4);
%% plot
figure,
plot(t,x(:,1)); xlabel('Time/s'),ylabel('Dimensionless displacement of primary structure')
figure,
plot(t,(x(:,2)); xlabel('Time/s'),ylabel('Dimensionless displacement of TMD ')
sys = ss(A,B,E,D);
figure,
bodeplot(sys(1,1))
bp.FrequencyScale = "linear";
[mag,phase,wout] = bode(sys(1,1));
mag = squeeze(mag);
phase = squeeze(phase);
fout = wout/(2*pi);
BodeTable = table(fout,mag,phase);
function dXdt = QZSdamper(t,A,B,X)
global alpha gamma mu Omega_1 Omega_2
w = 0.180;
F1 = 0.001/(Omega_1^2*L_0)*sin(w*t); %harmonic excitation
F_o = 2*alpha*(sqrt(1-gamma^2)+X(2)/L_0)*(sqrt(1/((X(2)/L_0)^2+2*sqrt(1-gamma^2)*X(2)/L_0+1))-1);
F = [F;F_o];
dXdt = A*X+B*(F);
end
5 Comments
Torsten
on 7 Oct 2024
Edited: Torsten
on 7 Oct 2024
I'm not a physician - thus I don't know how to non-dimensionalize time in a suitable way for your application. But as long as any quantity in your system of differential equations has a physical unit, the system is not dimensionless. And your derivatives are with respect to physical time t, not some kind of dimensionless time t', I guess.
In the definition of new parameters you probably mean
omega1 = sqrt(k1/m1), omega2 = sqrt(k2/m2)
instead of
omega1 = sqrt(k1^2/m1^2), omega2 = sqrt(k2^2/m2^2)
Then each term of equations (5) is dimensionless (except for the independent variable t).
Answers (2)
William Rose
on 7 Oct 2024
@Yu,
You are smart to add Sam Chak's name. He is a great and smart source of analysis and help.
There is more than one way you can create nondimensional coordinates for this system. It helps me to list all the physical constants in the systrem and relate them.
Before I do that, I suggest that the next time you post, you invest more time to provide a better diagram. This will make it easier for others to help you. The diagram on the left should be to the right of the diagram on the right, and it should be rotated 90 degrees, I think. The positive direction of x points up in the left diagram. The x direction is not shown on the right diagram, but I suspect it points to the right. Mass m2 and damper c2 appear in the discussion and in the equations, but neither m2 nor c2 appear in either diagram. And so on.
Is there a sign error in the equation for nonlinear force? You have
I think it should be . A better diagram would make it easier to know if the signs are correct.
Now let's list the physical constants of your system: masses m1, m2; lengths L0, a (or (a, h) or (h, L0)); dampers c1, c2; spring constants k1, kv, ko. You could write in terms of normalized length and normalized time. For length, you could normalize by L0. For time, you could choose to make the natural radian frequency, , of small oscillations of m1 about its rest point, equal to unity. To do this, you must compute in standard time units: , where knet is the sum of all the effectve spring constants acting on m1 at rest. It looks to me like , where is the angle of the ko springs at rest, and . I think the factor sine squared makes sense, but maybe there's an error or oversight in my mind's eye. Check my work.
4 Comments
Sam Chak
on 8 Oct 2024
Hi @Yu
The modeling process cannot be rushed, and the derived model must be verified. I advise young researchers to focus on "get things right before getting things done", as this can help avoid the need for repeated corrections. Therefore, before attempting to make the system dynamics dimensionless, please ensure that the model is correctly derived. You might also consider applying Lagrangian mechanics to see if you obtain the same mathematical model.
William Rose
on 8 Oct 2024
@Yu, The two diagrams in your recent comment are helpful.
You say in your recent comment "I define x2 as the relative displacement of m2 with respect to m1 to simplify the equations, as shown in my equations." I think this is a confusing definition, because the oblique springs are anchored at one end to the external world, not to m1. Therefore the force of the oblique springs (fo) depends on the position of m2 relative to the external world, and not on the4 position of m2 relative to m1. If you really want to define x2 as displacement of m2 relative to m1, then you need to revise your equations for f0 to reflect the fact that fo depends on (x1+x2) rather than on x2 only. You also need to revise your equation for , for the same reason. (Your equation for uses "x". I assume you mean x2. Please do not use "x" in your equations. Use only x1 and x2, for clarity.) For the reasons just listed, I encourage you to define x2 as position of m2 relative to the external world. I assume below that this is how x2 is defined.
Equation 3 in your original post is not correct for the system illustrated in the diagrams in your recent comment. The correct coupled equations for x1 and x2 are
10 Comments
David Goodmanson
on 11 Oct 2024
Edited: David Goodmanson
on 12 Oct 2024
Hello Yu,
In going dimensionless, in all cases I divided both sides of an equation by a constant, the same way as you are doing. In my comment (edited to number the equations), to go from (1) to (2) you plug in
d/dt = w0 d/du
(u being dimensionless) everywhere in (1), then divide both sides of the resulting equation by (m0 w0^2). The result is (2), in which the leading term is (M/m0).
In that term (M/m0) you you don't have to worry about factors of w0 any more because division by (m0 w0^2) takes care of it.
To get from (2) to (5), after substituting (4.2) into to (2), you have
M/m0) (d^2/du^2) x + C/(m0 w0) (d/du) x + K/(m0 w0^2) x =
= (1/(m0 w0^2)) 2 k0new h (1+x/h) (sqrt(h^2+a^2)/sqrt((h+x)^2+a^2) -1)
If you put in x = hy (which sort of implicitly assumes L0 = h) and divide both sides of the resulting equation by h, you get (5), with all terms dimensionless. The right hand side of (5) is
2 (k0new/(m0 w0^2)) (1+y) (sqrt(h^2+a^2)/sqrt((h^2(1+y)^2+a^2) -1)
You are of course right that there are still dimensions in the sqrt. But there is a ratio of two square roots and that ratio is dimensinless. I see now that you are interested in having dimensionless parameters everywhere as much as possible.. So in
sqrt(h^2+a^2)/sqrt((h^2(1+y)^2+a^2)
you can divide top and bottom by h, which is the same as dividing by h^2 inside the sqrt, to arrive at
sqrt(1+(a/h)^2)/sqrt(((1+y)^2+(a/h)^2)
where now the result is expressed in terms of a new dimensionless parameter (a/h).
I not commenting on whether (4.1) [your original nonlear force expression] is correct in tems of the system you are describing. It is just a good example of the process.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!