Can someone help me to obtain the FRF curve and plot the curve?

%%
% clear memory
clear all
clc

%% % Define material properties and beam dimensions % % %
%xx----------------------------------------------------------------------------------------xx%
width = 0.3;        % beam width
t1 = 0.1;           % beam thickness of first layer
t2 = 0.1;           % beam thickness of second layer
L = 1;              % beam length
poisson1 = 0.33;    % Poisson's ratio first layer
poisson2 = 0.33;    % Poisson's ratio second layer 
rho1 = 1700;        % density of material 1
rho2 = 1240;        % density of material 2
E1 = 2.2e11;        % modulus of elasticity for material 1
E2 = 7.9e19;        % modulus of elasticity for material 1
E3=7.9e9;           % modulus of elasticity for material 2

% % % beam load parameters % % %
%xx----------------------------------------------------------------------------------------xx%
n = 0;              % material gradation index
P = -100;           % uniform distributed load
Mb = 24.4;          % ratio of Young's moduli between two materials
E2 = E1 / Mb;

% % % Derived FGM parameters % % %
%xx----------------------------------------------------------------------------------------xx%
[C1, I01, I1, A01] = effectstifmass(width, t1, rho1, rho2, E1, E2, Mb, n, poisson1);

Unrecognized function or variable 'effectstifmass'.

[C2, I02, I2, A02] = effectstifmass(width, t2, rho1, rho2, E3, E3, Mb, n, poisson2);

% Mass parameters for FGM
numberElements = 40;
nodeCoordinates = linspace(0, L, numberElements + 1);
xx = nodeCoordinates'; 
x = xx';
elementNodes = zeros(size(nodeCoordinates, 2) - 1, 2);
for i = 1:size(nodeCoordinates, 2) - 1
    elementNodes(i, 1) = i;
    elementNodes(i, 2) = i + 1;
end
numberNodes = size(xx, 1);
GDof = 2 * numberNodes;

% % % Static analysis % % %
%xx----------------------------------------------------------------------------------------xx%
% Computation of the system stiffness, force, mass
[stiffness1, force, mass1] = formStiffnessMassTimoshenkoBeam(GDof, numberElements, ...
    elementNodes, numberNodes, xx, C1, P, rho1, I01, t1);
[stiffness2, ~, mass2] = formStiffnessMassTimoshenkoBeam(GDof, numberElements, ...
    elementNodes, numberNodes, xx, C2, P, rho1, I02, t2);

% % % Boundary conditions % % %
%xx----------------------------------------------------------------------------------------xx%
fixedNodeW = [1, 1];
fixedNodeTX = [1, 2];
prescribedDof = [fixedNodeW; fixedNodeTX + numberNodes];
stiffness = stiffness1 + stiffness2 + stiffness1;
mass = mass1 + mass2 + mass1;
K = stiffness;
m = mass;

% % % Free vibration analysis % % %
%xx----------------------------------------------------------------------------------------xx%
[modes, eigenvalues] = eigenvalue(GDof, prescribedDof, stiffness, mass, 0);
omega = sqrt(eigenvalues) * L * L * sqrt(rho1 * A01 / E1 / I01);
disp('First 3 dimensionless frequencies:')
disp(omega(1:3))
f = sqrt(eigenvalues);
disp('First 3 natural frequencies (Hz):')
disp(f(1:3))

% Drawing mesh and deformed shape
modeNumber = 3;
V1 = modes(:, 1:modeNumber);

% Normalize the mode shapes for better visualization
for i = 1:modeNumber
    V1(:,i) = V1(:,i) / max(abs(V1(:,i)));
end

% Plot the first 3 mode shapes
figure;
for i = 1:modeNumber
    subplot(3,1,i);
    plot(x, V1(1:numberNodes, i), 'LineWidth', 2);
    title(['Mode Shape ', num2str(i)]);
    xlabel('Beam Length (x)');
    ylabel('Normalized Displacement');
    grid on;
end

% Solve for static displacements
displacements = solution(GDof, prescribedDof, stiffness, force);
U = displacements;

% Max displacement
disp('Max normal displacement:')
disp(min(U(1:numberNodes)))

% % % Forced vibration analysis % % %
%xx----------------------------------------------------------------------------------------xx%
c1 = 1.5; 
c2 = 0.9;                       % damping coefficients
Ca = c1 * mass + c2 * stiffness;
F = @(t) force * sin(2* t);     % time-dependent load
D0 = zeros(GDof,1);             % initial displacement
V0 = zeros(GDof,1);             % initial velocity
dt = 0.005; T = 100; tt = 0:dt:T;
[D, V, A] = Newmark(mass, Ca, stiffness, F, D0, V0, dt, T);

% Plot forced vibration displacement at a node
figure;
plot(tt, D(21,:), 'b', 'LineWidth', 2);
title('Forced Vibration Displacement at Node 21');
xlabel('Time (s)');
ylabel('Displacement');
grid on;

% % % Export results to Excel % % %
%xx----------------------------------------------------------------------------------------xx%
xlswrite('mode vibration result.xlsx', [x(:), V1(1:numberNodes,:)]); 
xlswrite('normal displacement.xlsx', [x(:), U(1:numberNodes,:)]); 
xlswrite('forced displacement.xlsx', [tt', D(21,:)']);              

this is my current code

[Mathieu NOE]

hello

there are some functions in your code that you did not provide...

[Midhun]

can you please correct it. i am a beginner , i need to do the static and dynamic analyis of a layerd beam. could you please update me the code. also i dont konw to find the FRF

[dpb]

The source from which you obtained the above code would be the place to go first.

The Answers forum is not a general consulting platform; it is for specific questions about the use of MATLAB itself, it's the user's task to learn the basics of the problem to be solved.

[Mathieu NOE]

it's quite a bit complicated code for a beginner ... but we can help you

those lines are calling functions that are not part of the standard matlab package , so you must provided them to the forum (use the paper clip button) otherwise no one is able to run the code and extend it .

please provide the functions that appear in those lines :

[C1, I01, I1, A01] = effectstifmass(width, t1, rho1, rho2, E1, E2, Mb, n, poisson1);

[stiffness1, force, mass1] = formStiffnessMassTimoshenkoBeam(GDof, numberElements, ...
    
[modes, eigenvalues] = eigenvalue(GDof, prescribedDof, stiffness, mass, 0);
    
displacements = solution(GDof, prescribedDof, stiffness, force);

[D, V, A] = Newmark(mass, Ca, stiffness, F, D0, V0, dt, T);

[Sam Chak]

Hi @Midhun,

You need to be transparent with us. Since you posted the code without providing any background or explaining what you are trying to achieve with some of the functions, it is generally assumed that you wrote the code but encountered errors while running it.

If you inherited the original code from someone (or generated by Copilot) and modified it to suit your application, please specify which lines you understand and which lines you do not. The xlswrite() command is no longer recommended in current practice.

Please be specific and direct us to the lines in the code that perform the static and dynamic analyses. For example, does the following line perform the dynamic analysis?

[modes, eigenvalues] = eigenvalue(GDof, prescribedDof, stiffness, mass, 0);

Does this line perform the static analysis?

displacements = solution(GDof, prescribedDof, stiffness, force);

[Midhun]

can you please help me finf the FRF , like the command for obtaning the FRF and ploting the graph

[Sam Chak]

@Midhun, From your reply, I gather that you may be unable to provide the 5 functions that Mathieu requested to help you. In that case, you need to provide the mathematical formulas. Do you have the textbook, "Engineering Mechanics: Statics & Dynamics" by Hibbeler? We can likely guide you on how to implement those in MATLAB.

[dpb]

@Midhun downloaded or was given the unamed posted script but a search did not uncover the source for any of the missing functions nor do we know its provenance so have no way to find from whence it came and @Midhun isn't being forthright in volunteering any information nor shown any interest in providing any input from his end.

In line with blindly applying code without understanding what it is doing, see <Frequency Response Functions for Modal Analysis> in Signal Processing Toolbox.

Answer 1

[Sam Chak]

The beam data and load parameters are provided. Logically, by using the effectstifmass() function, the beam input data are used to obtain [C1, I01, I1, A01], which are required by the formStiffnessMassTimoshenkoBeam() function to return the mass, stiffness, and force matrices as outputs. These matrices typically constitute the MIMO system represented by the equation, , which can be expressed in the form of a state-space system. Subsequently, follow the guidance provided by @dpb to obtain the FRF for modal analysis.

Code snippet:

A   = [0  1;
      -1 -0.1];
B   = [0;
       1];
C   = [1, 0];
D   = 0*C*B;

Fs  = 1;
dt  = 1/Fs;
N   = 3000;
t   = dt*(0:N-1);
b   = 0.01;

u   = randn(1, N)/2;
u(2001:end) = 0;

y = 0;
x = [0;0];
for k = 1:N
    y(k) = C*x + D*u(k);
    x = A*x + B*u(k);
end

figure
plot(t, y), grid on
xlabel('Time')
ylabel('Displacement')

wind = hann(N/2);

%% Method 1
[frf, f]= modalfrf(u', y', Fs, wind, Sensor="dis");

%% Method 2
[b, a]  = ss2tf(A, B, C, D);
[ztf,fz]= freqz(b, a, 2048, Fs);

%% Plot
figure
hold on
plot(f, mag2db(abs(frf)))
plot(fz*Fs, mag2db(abs(ztf))), grid on
hold off
ylim([-20 60])
legend('Method 1', 'Method 2')
xlabel('Frequency')
ylabel('Magnitude')

[dpb]

@Sam Chak, you had more perserverence to actually read enough details in the beginning to see how much was/wasn't there.

I took the liberty to turn that effort into an Answer...hopefully @Midhun will express gratitude for others' efforts on his behalf.