# 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")