how to solve a 4 DOF problem with data input force

hy everyone,

I have troubles with a 4 DOF problem. In the specific, I have a 4 mass-spring problem with a damper and a force applied on the first mass, which values are available from an Excel sheet. I have to solve the differential equation for that system, but I don't know how to because the specific function of the input force is not defined and is a vector too. Could you help me? Also, I would like to know if it is possible to implement that problem in Simulink in order to check that the results are correct.

Best regards.

close all;clear all;clc

%DATA PROBLEM

A = xlsread('Database El Centro 1940_Primi 6 sec.xlsx');

time = A(:,2); %[s]

force = A(:,2); %[m/s^2]

E = 3*10^10; %[Pa]

I = (3.1)*10^(-3); %[m^4]

L = 3; %[m]

m = (8.04)*10^3; %[kg]

k = 48*E*I/(L^3); %[N/m]

csi = 0.2;

c = 2*csi*sqrt(m*k); %[Ns/m]

%MASS & STIFFNESS MATRICES

M = [m 0 0 0;0 m 0 0;0 0 m 0;0 0 0 m];

C = [c 0 0 0;0 0 0 0;0 0 0 0;0 0 0 0];

K = [2*k -k 0 0;-k 2*k -k 0;0 -k 2*k -k;0 0 -k k];

f = zeros(length(M),length(time));

f(1) = force';

dx2_dt2 = @(t,x)[[x(5),x(6),x(7),x(8)]'; -C*[x(5),x(6),x(7),x(8)]' ...

- K*[x(1),x(2),x(3),x(4)]' + f];

odeopt = odeset('RelTol',0.001,'AbsTol',0.001,'InitialStep',0.02,...

'MaxStep',0.02);

x_0 = zeros(1,length(M));

x_0_dot = zeros(1,length(M));

[time,x] = ode45(dx2_dt2,[0 6],[x_0 x_0_dot],odeopt);

Answer 1

[Sam Chak]

Hi @LUCA PATRIARCHI

Since the earthquake force is externally time-dependent (different step size from the solver), I have created a function named LucaP that accepts four input arguments: t, x, ft, and f, so that the force can be interpolated inside the ODE function.

function dxdt = LucaP(t, x, ft, f)

dxdt = zeros(8,1);

E = 3*10^10; %[Pa]

I = (3.1)*10^(-3); %[m^4]

L = 3; %[m]

m = (8.04)*10^3; %[kg]

k = 48*E*I/(L^3); %[N/m]

csi = 0.2;

c = 2*csi*sqrt(m*k); %[Ns/m]

%MASS & STIFFNESS MATRICES

M = [m 0 0 0;0 m 0 0;0 0 m 0;0 0 0 m];

C = [c 0 0 0;0 0 0 0;0 0 0 0;0 0 0 0];

K = [2*k -k 0 0;-k 2*k -k 0;0 -k 2*k -k;0 0 -k k];

f = interp1(ft, f, t); % Interpolate the data set (ft, f) at times t

dxdt(1:4) = [x(5),x(6),x(7),x(8)]';

dxdt(5:8) = -C*[x(5),x(6),x(7),x(8)]' - K*[x(1),x(2),x(3),x(4)]' + f;

end

To solve the ODE using ode45, you need to specify the function handle so that it passes the vectored values for ft and f to LucaP for the interpolation job. Note that in the original script, there was a typo on time = A(:,2);.

A       = xlsread('Database El Centro 1940_Primi 6 sec.xlsx');
ft      = A(:,1);           % time vector extracted from Excel
f       = A(:,2);           % force vector extracted from Excel
tspan   = [0 6];            % simulation time span
x_0     = zeros(1,4);
x_0_dot = zeros(1,4);
ic      = [x_0 x_0_dot];    % initial condition

odeopt  = odeset('RelTol', 0.001, 'AbsTol', 0.001, 'InitialStep', 0.02, 'MaxStep', 0.02);
[t, x]  = ode45(@(t, x) LucaP(t, x, ft, f), tspan, ic);

Result:

Note: The four velocity states are toggled off because the signals are highly oscillatory.

[LUCA PATRIARCHI]

Hy @Sam Chak,

thank you for your answer, I really appreciated that. For what concerns a Simulink simulation of the same problem, could you help me? I am asking it only because I believe that having a controller for future problems would be usefull.

Best regards.

[Sam Chak]

Hi @LUCA PATRIARCHI

Don't mention it. You can post a new question if you have troubles in Simulink. Since your model is a linear system, there are relatively many control design tools you can apply. More importantly, you should find more details related to Tuned Mass Dampers or Seismic Dampers. All the best in your future endeavors.

https://www.sciencedirect.com/search?qs=Tuned%20mass%20damper