Transient Beam Model

Build a transient model and compute the impulse response.

First construct the beam and specify the Young's modulus, Poisson's ratio, and mass density of steel. Specify the tip of the beam using the addVertex property.

gm = multicuboid(0.1,0.005,0.005);
E = 210E9;
nu = 0.3;
rho = 7800;
TipVertex  = gm.addVertex('Coordinates',[0.05,0,0.005]);
firstNF = 2639;
Tfundamental = 2*pi/firstNF;

Use createpde (Partial Differential Equation Toolbox) to construct the transient model and generate the mesh using the generateMesh (Partial Differential Equation Toolbox) command.

transientModel = createpde('structural','transient-solid');
transientModel.Geometry = gm;
msh = generateMesh(transientModel);

Assign structural properties for the steel beam with the structuralProperties (Partial Differential Equation Toolbox) command and fix one end using structuralBC (Partial Differential Equation Toolbox).

structuralProperties(transientModel,'YoungsModulus',E,'PoissonsRatio',nu,'MassDensity',rho);
structuralBC(transientModel,'Face',5,'Constraint','fixed');

Visualize the beam geometry using pdegplot (Partial Differential Equation Toolbox).

pdegplot(transientModel,'VertexLabels','on');
title('Beam model')

Apply force on beam tip for 2% of fundamental period of oscillation (impulse) using structuralBoundaryLoad (Partial Differential Equation Toolbox).

Te = 0.02*Tfundamental;
structuralBoundaryLoad(transientModel,'Vertex',TipVertex,...
   'Force',[0;0;-100],'EndTime',Te);

Set initial conditions for the transient beam model using structuralIC (Partial Differential Equation Toolbox).

structuralIC(transientModel,'Velocity',[0;0;0]);

Compute the impulse response by solving the transient beam model.

ncycles = 10;
tsim = linspace(0,ncycles*Tfundamental,30*ncycles);
R1 = solve(transientModel,tsim);

Visualize the oscillations at the tip of the beam.

plot(tsim,R1.Displacement.uz(TipVertex,:))
title('Vertical displacement of beam tip')
legend('Transient model')
xlabel('Time')
ylabel('Displacement')

The next step is to extract the sparse matrices from the finite element beam model.

Sparse Model from the Transient Model

The first step to extract the M,K,B,F matrices from the finite element beam model to modulate the force to compute the B matrix. Use assembleFEMatrices (Partial Differential Equation Toolbox) to extract the sparse matrices for the second-order beam model.

F1 = -100;
structuralBoundaryLoad(transientModel,'Vertex',TipVertex,'Force',[0;0;F1]);
state.time = 0;
data1 = assembleFEMatrices(transientModel,'nullspace',state);
F2 = 100;
structuralBoundaryLoad(transientModel,'Vertex',TipVertex,'Force',[0;0;F2]);
data2 = assembleFEMatrices(transientModel,'nullspace',state);

The resultant model is defined by the following set of equations:

M = data1.M;
K = data1.Kc;

The F matrix consists of x, y and z displacements at the TipVertex node.

N = size(msh.Nodes,2);
F = data1.B;                  
F = F(TipVertex+[0 N 2*N],:);

Map the B matrix relating the force to the input.

B = (data2.Fc-data1.Fc)/(F2-F1);

Obtain the sparse model of the beam using mechss.

sys = mechss(M,[],K,B,F(3,:));

Use spy to visualize the sparsity of the mechss model sys.

spy(sys)

Linear Analysis

Enable parallel computing and choose 'tfbdf3' as the DAE solver.

sys.SolverOptions.UseParallel = true;
sys.SolverOptions.DAESolver = 'trbdf3';

Use bode to compute the frequency response of this model.

w = logspace(2,6,1000);
bode(sys,w), grid

using parallel (mechss)

title('Frequency response from force to tip vertical displacement')

Next use lsim to compare the impulse response with the approximations obtained from the transient beam model. To limit error due to linear interpolation of force between samples, use a step size of Te/5. The force applied at the beam tip is $0<=t<=\mathrm{Te}$.

h = Te/5;  
t = 0:h:ncycles*Tfundamental;
u = zeros(size(t));    
u(t<=Te) = -100;

y = lsim(sys,u,t);
plot(t,y,tsim,R1.Displacement.uz(TipVertex,:))
title('Comparison of full-model simulated response with transient model approximations')
legend('Simulated model','Transient model')
xlabel('Time')
ylabel('Displacement')

The transient model approximations match well with the full-model simulation using the TRBDF3 DAE solver.