I would like to validate some code that I have that seems to be unsuccessful.

A uniform beam of length (L), Elastic Modulus (E), Second Moment of Area (I), mass per unit length (rho), and damping ratio (zeta) is subjected to the arbitrary time-varying forces fe(t) at x = xe, fa(t) at x = xa, and time-varying moment/torque ma(t) at x = xa. Euler-Bernoulli assumptions are made and the Ritz-Galerkin modal approximation (assumed-modes method) is used to describe the dynamic motion of the beam.

For a given L, E, I, rho, zeta, fe(t), fa(t), ma(t), xe, xa...
I want to model the deflection of the beam y(x,t) with both spatial (x) and temporal (t) dimensions taken into account.

Does anyone have some code that could create this model for me? I am including a simple version of the code I am currently using below. However, myself and the people I am working with are a little bit too deep into the problem that we are using the code for and haven't been able to find any mistakes. Therefore, it might be more helpful to ignore my code and start from scratch to have an unbiased approach. Then compare the output figures from my code and yours.

```
clear;
nmodes = 5; % Number of modes used in the modal approx.
L = 1; % Beam length
E = 1; % Beam Elastic Modulus
Iner = 1; % Beam moment of inertia
rho = 1; % Beam mass per unit length
zeta = 0.001; % Beam damping ratio
xe = 0.75*L; % Location of force fe(t), normalized to beam length
xa = 0.2*L; % Location of force fa(t) and moment ma(t), normalized to beam length
[phixe,~,PHI,N,S] = ex_modeshape(nmodes,xe,L,rho,E,Iner);
[phixa,phixa_x,~,~,~] = ex_modeshape(nmodes,xa,L,rho,E,Iner);
C = 2*zeta*sqrt(N.*S); % Modal damping
% System Equation of Motion: AQ = F
% A (nmodes x nmodes): characteristic matrix
% Q (nmodes x 1): modal coordinates in Laplace domain... Y(x,s) = sum(Q(s).*PHI(x))
% F (nmodes x 1): generalized force vector
syms s
for ii = 1:nmodes
A(ii,ii) = N(ii)*s^2 + C(ii)*s + S(ii); % Note: Laplace domain
end
%% To find the steady-state amplitudes of the beam under HARMONIC EXCITATION
% fe(t) = Fe*sin(om*t)
% fa(t) = Fa*sin(om*t)
% ma(t) = Ma*sin(om*t)
Fe = 1;
Fa = 2;
Ma = 0.5;
om = 10; % Excitation frequency
s = om*i;
F =Fe*phixe + Fa*phixa + Ma*phixa_x;
Q = eval(A\F); % CAUTION: This becomes singular if the system is harmonically excited at any of its natural frequencies
syms x real
Yi(x) = sum(Q.*PHI(x)); % y(x,s=w*i). Note: still in frequency domain
XSIM = 0:0.01*L:L; % Specify density of points along beam length to plot
for ll = 1:length(XSIM)
xsim = XSIM(ll);
Y(ll) = Yi(xsim);
Y_amp(ll) = abs(Y(ll))*sign(real(Y(ll))); % y(xsim,s=w*i)
end
figure();
plot(XSIM,Y_amp);
grid on;
xlabel('x/L');ylabel('|Y|_{ss}');
title('Steady State Amplitudes');
```

ex_modeshape.m function:

```
function [phix,phix_x,PHI,N,S] = ex_modeshape(nmodes,xeval,L,rho,E,Iner)
% ex_modeshape.m is a function that outputs...
%
% phix (nmodes x 1): Mode shapes for fixed-free beam evaluated @ xeval
%
% phix_x (nmodes x 1): Partial derivatives of fixed-free mode shape
% functions evaluated @ xeval
%
% PHI (nmodes x 1 symfun): A vector of the actual mode shape functions
%
% N (nmodes x 1): The modal masses, under mass-normalization
%
% S (nmodes x 1): The modal stiffnesses, under mass-normalization
B = [1.875104;...
4.694091;...
7.85475744;...
10.995540735;...
14.13716839105;...
17.278759532088;...
20.42035225104125]/L; % The modal constants for the first 7 modes
syms x real
for ii = 1:nmodes
if ii > 7
B(ii) = (2*ii-1)*pi/(2*L); % After the first 7 modes, the approximation for the modal constants becomes sufficiently precise.
end
phi(x) = cos(B(ii)*x) - cosh(B(ii)*x) + ((sin(B(ii)*L) - sinh(B(ii)*L))/(cos(B(ii)*L) + cosh(B(ii)*L)))*(sin(B(ii)*x) - sinh(B(ii)*x)); % The actual mode shape function
b(ii,1) = sqrt(eval(1/(rho*int((phi(x))^2,x,0,L)))); % Scaling factor calculated under mass-normalization
fee(ii,1) = b(ii,1)*phi(x); % The scaled (mass-normalized) mode shape
phi_x(x) = b(ii)*diff(phi(x),x,1); % d_phi/d_x
N(ii,1) = 1; % Modal mass
S(ii,1) = B(ii)^4*E*Iner/rho; % Modal stiffness equals natural frequency squared under mass-normalization
phix(ii,1) = b(ii)*eval(phi(xeval));
phix_x(ii,1) = eval(phi_x(xeval));
end
PHI(x) = fee;
end
```

This is for a simplified case with harmonic external forces and moment, all in phase with one another. Ideally, these external excitations are more generalized, or at least include some phase shift between them.