Modal Superposition Method for Structural Dynamics Problem

This example shows how to solve a structural dynamics problem by using modal analysis results. Solve for the transient response at the center of a 3-D beam under a harmonic load on one of its corners. Compare the direct integration results with the results obtained by modal superposition.

Modal Analysis

Create the geometry and plot it with the edge and vertex labels.

gm = multicuboid(0.05,0.003,0.003);
pdegplot(gm,EdgeLabels="on", ...
            VertexLabels="on", ...
            FaceAlpha=0.3);
view([95 5])

Create an femodel for modal structural analysis and include the geometry.

model = femodel(AnalysisType="structuralModal", ...
                 Geometry=gm);

Generate a mesh.

model = generateMesh(model);

Specify Young's modulus, Poisson's ratio, and the mass density of the material.

model.MaterialProperties = ...
    materialProperties(YoungsModulus=210E9, ...
                       PoissonsRatio=0.3, ...
                       MassDensity=7800);

Specify minimal constraints on one end of the beam to prevent rigid body modes. For example, specify that edge 4 and vertex 7 are fixed boundaries.

model.EdgeBC(4) = edgeBC(Constraint="fixed");
model.VertexBC(7) = vertexBC(Constraint="fixed");

Solve the problem for the frequency range from 0 to 500000. The recommended approach is to use a value that is slightly smaller than the expected lowest frequency. Thus, use -0.1 instead of 0.

Rm = solve(model,FrequencyRange=[-0.1,500000]);

Transient Analysis

Change the analysis type to transient structural.

model.AnalysisType = "StructuralTransient";

Define a sinusoidal load function, sinusoidalLoad, to model a harmonic load. This function accepts the load magnitude (amplitude), the location and state structure arrays, frequency, and phase. Because the function depends on time, it must return a matrix of NaN of the correct size when state.time is NaN. Solvers check whether a problem is nonlinear or time-dependent by passing NaN state values and looking for returned NaN values.

function Tn = sinusoidalLoad(load,location,state,Frequency,Phase)
if isnan(state.time)
    Tn = NaN*[location.nx location.ny location.nz];
    return
end
if isa(load,"function_handle")
    load = load(location,state);
else
    load = load(:);
end
% Transient model excited with harmonic load
Tn = load.*sin(Frequency.*state.time + Phase);
end

Apply a sinusoidal force on the corner opposite the constrained edge and vertex.

Force=[0 0 10];
Frequency = 7600;
Phase = 0;
forcePulse = ...
    @(location,state) sinusoidalLoad(Force, ...
                                     location,state, ...
                                     Frequency,Phase);
model.VertexLoad(5) = vertexLoad(Force=forcePulse);

Specify the zero initial displacement and velocity.

model.CellIC = cellIC(Velocity=[0;0;0],Displacement=[0;0;0]);

Specify the relative and absolute tolerances for the solver.

model.SolverOptions.RelativeTolerance = 1E-5;
model.SolverOptions.AbsoluteTolerance = 1E-9;

Solve the model using the default direct integration method.

tlist = linspace(0,0.004,120);
Rd = solve(model,tlist);

Now, solve the model using the modal results.

tlist = linspace(0,0.004,120);
Rdm = solve(model,tlist,ModalResults=Rm);

Interpolate the displacement at the center of the beam.

intrpUd = interpolateDisplacement(Rd,0,0,0.0015);
intrpUdm = interpolateDisplacement(Rdm,0,0,0.0015);

Compare the direct integration results with the results obtained by modal superposition.

plot(Rd.SolutionTimes,intrpUd.uz,"bo")
hold on
plot(Rdm.SolutionTimes,intrpUdm.uz,"rx")
grid on
legend("Direct integration", "Modal superposition")
xlabel("Time");
ylabel("Center of beam displacement")