Nondimentionalize the 2 DOF equations

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

[Yu]

@Sam Chak

[Torsten]

I think mu = m1/m2 is incorrect. It should be mu = m2/m1.

And in system (5) each term has unit s^2 - thus your system is not dimensionless. You must write your equation not in x1(t) and x2(t), but in x1(t'*T)/L0 and x2(t'*T)/L0 where T is some characteristic time (e.g. the endtime) (0<=t'<=1 dimensionless, defined as t/T).

And the easiest way to compare whether you did everyting right is to compare the results from (3) and (5). They should agree.

[Yu]

@Torsten Sorry for the mistakes. I have modified the equations just now.

Thank you for you suggestions. But I still have two questions:

1. Do I just need to divide the T (e.g. endtime) for achieving completely dimensionless?
2. I have divided the ω (/s^2) in the equations. Is this not used to eliminate the unit s^2-?

[Torsten]

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).

[Yu]

Hi @Torsten, I am sorry they are indeed mistakes :( They should be

and

In this case, equation (5) is dimensionless (except for the independent variable t) as you mentioned.

Answer 1

[William Rose]

@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.

[Yu]

Hi William,

I really appreciate your valuable comments. As it has been night in my place, I would like to read your comments carefully tomorrow. I just followed the suggestions from Sam Chak in my previous post. As there is a messy description of my questions.

I will make the diagram better. Thank you!

I apologize for the missing information. It is positive when the moves to right.

I did forget to post the equations derivation in more detailed way.

As you can see, the three-spring system in the left diagram will provide nonlinear stiffness force. In my assumption, the oblique srpings () is pre-compressed. When m_2 moves to the right, the force would be . Therefore, the direction of oblique stiffness force would be opposite to the direction of vertical stiffness force. As for my equations, it is basically in the form of . Thus kx term should be positive when assuming the vertical stiffness force > oblique stiffness force (it is supposed to be).

The f_non has been verified:

I am not sure if I clarify in a clear way. If you are confused about my description, please let me know.

Best wishes,

Yu

[William Rose]

@Yu,

Thank you for your response and for the diagram. You refer above to "vertical" and "oblique" forces, and you say "the direction of oblique stiffness force would be opposite to the direction of vertical stiffness force". This suggests that "oblique" and "vertical" are opposite directions. Please explain forces and directions relative to the +x axis. I realize that forces orthogonal to the x-axis cancel out, since the oblique springs are assumed to be identical but opposite.

You said you would provide a better diagram, but you have not. I am trying to imagine the physical arrangement that accounts for equation 3 in your original post. If we convert equation 3 from a matrix equation to a pair of coupled equations, we get:

Here are two diagrams I can imagine, and the associated differential equations. Neither diagram's differential equations match your equation 3. In my equations below, should be . Force acts on m2 and is due to oblique springs.

I assume the oblique springs Hookean, and that the nonlinear behavior arises soley because they are oblique. Therefore, in the limit as , the nonlinear springs become linear, because in this case they are aligned with the x-axis. Is that your understanding also?

In your comment, you include a plot of fnon versus xnon. What is xnon? Does it equal x2? Indicate xnon and x2 on your diagram, and explain if they are different. I understand (from your original posting) that, on the plot of fnon, .

In your original post, you say(your equation 1). Does "x" in equation 1 refer to x1 or x2? I assume from equation 1 that fnon is the force due to the network of kv plus two oblique springs. It appears from equation 3 that fnon acts on m2. But equation 1 indicates that when x (i.e. x2) becomes more positive, the kv component of fnon also becomes more positive. For a stable system, the force from spring kv should become more negative, or less positive, when x becomes more positive, so that spring kv provides a restoring force. This is why I think the sign for the component of fnon, in equation 1, is incorrect. The sign of the component is correct in equation 3.

For the plot of fnon versus xnon, please confirm that fnon includes both the kv component AND the ko components. It might be more instructive to plot the ko component of force only, versus x2, since kv is a simple linear spring which is well understood.

Equation 2 in your original post indicates that the force from oblique springs equals zero only when x=-h. The position x=-h corresponds to the spirings being oriented perpendicular to the direction of motion. In other words, equation 2 suggests that the oblique springs are unstressed, at their resting length, when x=-h. Your system diagram indicates that when x=0, the springs are longer than when x=-h. This means the springs are extended, not compressed, when x=0, which is contrary to what you said in your comment.

I assume that m1 and m2 are initially at equilibrium and that the masses are at positions x1=0 and x2=0 respectively. Do you agree? You said the oblique springs are compressed in the initial equlibrium state. Then k1 and maybe kv must also be compressed, for the system to be at equilibrium. If you think the initial compression of the oblique springs is important, then the equations for all springs may need to include the unstressed length of each, and you should confirm that the spring forces sum to zero, on each mass, when x1=0 and x2=0 - that is, if the system is supposed to be at equilibrium at x1=0 and x2=0.

I am not sure why you want to find non-dimensional equations for the system. If it were me, I would want to make sure I have the equations right, and solve the system, using the real physical units. Then I mihght try to gain deeper insight thorugh dimensional analysis.

[Yu]

Hi @William Rose,

I really appreciate your detailed explantion. I would admit that I made a fatal mistake. The system should satisfy the physical meaning. When x1=0 and x2=0, the system should have 0 external force which means the oblique springs should not be compressed in the initial equlibrium state.

However, I’m a bit confused about Equation 2 from my original post. When x2 (x) =0, the equation in the red rectangle below should also equal zero. Yet, if I assume there is compression at the initial position (x=0), it shouldn’t be zero.

I am sorry for the late response. I made two diagrams for explaining my thought. I hope they are better now.

In fig. 1, the three-spring system replaces the linear spring (if it could function that way). When m1 and m2 are at their initial position, there is no force acting on them, meaning the springs are neither extended nor compressed. I define x2 as the relative displacement of m2 with respect to m1 to simplify the equations, as shown in my equations. When a displacement occurs, induced by f(t), the force generated by the three-spring system is expected to be:

(1)

(2)

where the extented lenght and .

It is shown that is 0 when x_2=0 or x_2=-h. When x = 0, the system is at the initial position there is no displacement and no force. When x = -h, the oblique spring is horizontal. I think this would be more reasonable.

Additionally, I believe there is another configuration for the three-spring system, as shown in fig. 2. In the initial position, the oblique springs are horizontal, allowing them to be compressed to their length L0. When a displacement occurs, the restoring force can be expressed as:

where , L0 is the original length of the oblique spring and l is the compressed lenght of the oblique spring.

Regarding the plot of the three-spring system, it incorporates a vertical spring component. The objective of the system is to achieve zero stiffness, as illustrated in the plot. The stiffness force will be zero when x2 is at a specific position.

I aim to find the dimensionless form of the system because there are numerous parameters that need to be determined, such as h, a, and L0, which are all included in the stiffness force equations. Since the dimensional ratios are related to the stiffness ratios, I would like to understand how the stiffness ratio influences the primary structure. When I nondimensionalize the , I also need to nondimensionalize the whole equation.

I would try your suggestions to make sure the equations correct.

Thank you again for your detailed explanations. They are really helpful.

Best wishes,

Yu

[Sam Chak]

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.

Answer 2

[William Rose]

@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

[William Rose]

I am having trouble posting my replies. Please bear with me.

[William Rose]

Equation 3 in your original post is incorrect, in my opinion. The correct coupled equations for x1 and x2 are

...

which you can write in matrix form. Check my work.

You say in your recent comment "The objective of the system is to achieve zero stiffness, as illustrated in the plot. The stiffness force will be zero when x2 is at a specific position." Please define "zero stiffness" and "stiffness force". Does this objective apply to x1 or to x2 or to both? I suspect that what you mean is that ... when x2=0, and ... when x2=0. But I could be wrong.

You use L and {\ell} in your recent post. Please explain with equations what the difference is.

I will discuss in a later post. No time now.

[Yu]

Hello William,

Thank you for you comments. I would agree that when the oblique springs are connected to the external word, the stiffness force would depend on (x1+x2). Basically, the stiffness force is defined as kx (stiffness multiple the displacement), and stiffness is k.

The 'zero stiffness' is that as you can see in the plot below.

fnon is the forces including the vertical and oblique spring forces. Similarly, Knon (by differentiating the fnon) is the stiffness including the vertical and oblique spring stiffness. xnon is the displacement of the m2. This system can achieve the 'zero stiffness' when . This position for nondimensional displacement should be . (For better observation, I make the x-axis in below plot positive).

I am not sure if you mean and l in the fig.2 in my previous post. is the original length of the oblique spring in fig.2 and l is the compressed length of the oblique spring.

If the oblique spring are connected to the external world, the stiffness force from the oblique spring will be only applied to the m2.

I am sorry that it seems I can not see some contents you presented.

Also, I have double-checked my 2DOF equations without with the results from simulink. They have a good agreement (@Torsten). I suppose the nonlinear stiffness force would be the critical problem. I am also now further verifying my equations based on Lagrange equations.

[William Rose]

@Yu,

I am finding it impossible to post comments that contain equations written with the Matlab Answers equation editor. I tried again, still no luck. I also tried screenshot of the equations and inserting the image, but that too does not post.

[Yu]

@William Rose

I have on idea what happened. I used screenshot in my previous posts.

I am trying to compare the results between dimensionless form and the real system by following the suggestions from Torsten.

[David Goodmanson]

Hello Yu,

I agree with WR here. First that the equations have to be complete and correct before even thinking about dimensionality. Second, it makes sense for displacement x2 of m2 to be relative to the fixed lab frame rather than relative to m1. Then for the linear part of the problem the stiffness matrix K is no longer diagonal but it's still easy, being

[k1+k2, -k2
  -k2,   k2]

(similar for C) and the harder, nonlinear part of the force is easier to calculate.

After the equations are set up, including unstretched spring lengths, it will be useful to have the static equilibrium position, where the total potential energy for all the springs is at a minimum. That can be done numerically if necessary.

As to dimensionality, you have a system

M (d^2/dt^2) x  + C (d/dt) x + K x = F(t,x)               (1)

where x is the column vector [x1;x2] and F is the nonlinear force contribution due to the oblique springs. (In your case F is not an explicit function of t, but in general it could be). For dimensionless equations you need characteristic parameters m0,L0,t0. The characteristic mass can be one of the masses, say m0 = m1, and the characteristic time can be t0 = 1/w0 = sqrt(m0/k0), where k0 is a suitable choice from K, but that's it. Bringing in a parameter from C, say c0 = c1, provides a useful dimensionless constant, the damping ratio c0/(w0 m0) [ or c0/(2 w0 m0) as in your eqn (6) ], but no L0. For that you need the force term.

Introducing dimensionless time u, with t = u/w0, d/dt = w0 d/du gives

(M/m0) (d^2/du^2) x  + C/(m0 w0) (d/du) x + K/(m0 w0^2) x = F(t,x)/(m0 w0^2)    (2)
% all terms are proportional to length

There are a couple of ways to get a the characteristic length. If the force is defined by a specific value F0 as, for example, F = F0 cos(wt) then

L0 = F0/(m0 w0^2)                                                        (3)

That's not the case you have though. Instead there is a bunch of geometry and spring constants. However, the nonlinear force can be expressed as F = k L f(x2). Here k and L are a characteristic spring constant and distance for the nonlinear part, and f(x) is dimensionless function of x that varies within a few orders of magnitude around 1. Taking your orginal example you have

F = 2 k0new (h+x) (sqrt(h^2+a^2)/sqrt((h+x)^2+a^2) -1)                 (4.1)
% here k0new is not the same k0 as was used before
= 2 k0new h (1+x/h) (sqrt(h^2+a^2)/sqrt((h+x)^2+a^2) -1)               (4.2)

An obvious choice is L0 = h and after going to x = hy with a dimensionless y, the final result is the dimensionless

(M/m0) (d^2/du^2) y  + C/(m0 w0) (d/du) y + K/(m0 w0^2) y
        = 2 k0new/(m0 w0^2) (1+y) (sqrt(h^2+a^2)/sqrt((h^2(1+y)^2+a^2) -1)    (5)

[Yu]

Dear @David Goodmanson,

Many thanks for your comments.

I have modified my equations and compare the results (regarding the relative frame and fixed frame). The results shows that they have some difference in the amplitude.

Thank you for clearing things up. I still have some questions.

1. From your first equations to the final equations, I think they are really clear. But I am not sure if I am correct. I think the final equation is similar with mine.
2. And it seems that you were nondimensionlizing the terms in the equations separately. I nondimensionlized the equations by dividing 'term' on both side.

% Subsitute the x = hy
(M/m0)(d^2/du^2)(hy) % I think here should be w0 as t0=u/w0. But there should be also unit s- left
C/(m0 w0)(d/du)(hy) % I think it is similar with the damping term of my equations
K/(m0 w0^2)(hy) % This also make sense for me
F = 2k0new h/ (m0 w0^2)(1+y) (sqrt(h^2+a^2)/sqrt((h+x)^2+a^2) -1) % I am a bit confused with this part. 
% There are still some dimensions in the part of 'sqrt'. 
% h can be eliminated by divided h on both side?

I would clarify the reason why I decide to nondimensionlize. I use L0 to nondimensonlize the equations.

1. As you can see the nonlinear stiffness force equation, there are many geometric parameters. I would like to present the dynamic behaviour of the system. Therefore, I would present a more 'generic' system without any specific values but ratios. Like the dimensional ratios and stiffness ratio (k0/kv).
2. If you see the configuration of the nonlinear stiffness system, there are some relationships: , . Hence, if I use L0 as the characteristic length, the equation for the nonlinear force would be:

F = 2k0new /(m0 w0^2)(1+y)(1/sqrt(2 sqrt(1-gamma^2) y +1+ y) - 1) % y = x/L0

I am not sure whether can I choose L0 as the characteristic length. And why there is no w0 in the M term.

Best wishes,

Yu

[Yu]

Hi @David Goodmanson,

I review your third equation again. I am not sure how you get this form. It seems that you divided different terms for this equations.

[David Goodmanson]

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.

[Yu]

Hi @David Goodmanson,

Thank you for your patient explantionation. According to your comments, I reorganized my equations to achieve dimensionless. Unfortunately, I found some difference and checked the process but had no clue. Thereofore, I would present the results for some support. Below is the figs about the dimensional and dimensionless displacement of the primary structure. I just removed the external excitations and considered the free vibration with the initial conditions.

It is obvious that they have different frequency. But I can confirm that the parameters are same.

There are equations for these two model in matlab.

As you can see, I divided for these equations. Where , and ,, . Please note k_2 is the vertical spring stiffness of the three spring system. xi_1 and xi_2 = .

For the both cases, I set the parameters. like damping ratio and mass ratio to be same. But the results are really different.

Best wishes,

Yu