Define Parameters for Structural Analysis

Create a square cross-section beam geometry.

gm = multicuboid(0.05,0.003,0.003);

Plot the geometry, displaying face and edge labels.

figure
pdegplot(gm,FaceLabels="on",FaceAlpha=0.5)
view([71 4])

figure
pdegplot(gm,EdgeLabels="on",FaceAlpha=0.5)
view([71 4])

Add a vertex at the center of face 3.

centerVertex = addVertex(gm,Coordinates=[0.025 0.0 0.0015]);

Create an femodel object for transient structural analysis and include the geometry in the model.

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

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

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

Fix one end of the beam.

model.EdgeBC([2 8 11 12]) = edgeBC(Constraint="fixed");

Generate a mesh.

model = generateMesh(model);

Apply a sinusoidal concentrated force in the z-direction on the new vertex. First, define a sinusoidal load function, sinusoidalLoad, to model a harmonic load. This function accepts the load magnitude (amplitude), 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)
    normal = [location.nx location.ny];
    if isfield(location,"nz")
        normal = [normal location.nz];
    end
    Tn = NaN*normal;
    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

Now apply the force on the new vertex.

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

Specify zero initial conditions.

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

Reduce Model

Specify the fixed and loaded boundaries as structural superelement interfaces by creating a romInterface object for each superelement interface. The reduced-order model technique retains the DoFs on the superelement interfaces while condensing all other DoFs to a set of modal DoFs. For better performance, use the set of edges bounding face 5 instead of using the entire face.

romObj1 = romInterface(Edge=[2 8 11 12]);
romObj2 = romInterface(Vertex=centerVertex);

Assign a vector of interface objects to the ROMInterfaces property of the model.

model.ROMInterfaces = [romObj1,romObj2];

Reduce the structure, retaining all fixed interface modes up to 5e5.

rom = reduce(model,FrequencyRange=[-0.1,5e5]);

Simulate Transient Dynamics Using ROM

Next, use the reduced-order model to simulate the transient dynamics. Use the ode15s function directly to integrate the reduced system of ordinary differential equations. Take the loaded and modal DoFs for time-integration, and leave the fixed DoFs aside because the solution remains zero for those DoFs.

Working with the reduced model requires indexing into the reduced system matrices rom.K and rom.M. The arrangement of DoFs in reduced system is such that the physical DoFs corresponding to retained interfaces appear first followed by the generalized model DoFs. DoFs in a structural problem correspond to translational displacements. If the number of mesh points in a model is Nn, then the software assigns the IDs to the DoFs as follows: the first 1 to Nn are x-displacements, Nn+1 to 2*Nn are y-displacements, and 2Nn+1 to 3*Nn are z-displacements. Only the subset of these 3*Nn DoFs corresponding to ROMInterfaces is retained in the reduced model. The reduced model object rom contains these IDs for the retained DoFs in rom.RetainedDoF.

Create a function that returns DoF IDs given node IDs and the number of nodes.

getDoF = @(x,numNodes) [x(:); x(:) + numNodes; x(:) + 2*numNodes];

Find the node at the loaded vertex.

loadedNode = findNodes(rom.Mesh,"region",Vertex=centerVertex);

Find the DoF of the loaded nodes using the helper function getDoF.

numNodes = size(rom.Mesh.Nodes,2);
loadDoFs = getDoF(loadedNode,numNodes);

Knowing the DoF IDs for the given node IDs, use rom.RetainedDoF and the intersect function to find the required indices corresponding to those DoF in the reduced matrices.

[~,loadNodeROMIds] = intersect(rom.RetainedDoF,loadDoFs);

In the reduced matrices rom.K and rom.M, generalized modal DoFs appear after the retained DoFs. Find the indices of modal DoFs in rom matrices.

modelDoFIDs = ((numel(rom.RetainedDoF) + 1):size(rom.K,1))';

Find the indices for the ODE DoFs in reduced matrices. Because fixed-end DoFs are not a part of the ODE system, these indices are as follows.

odeDoFs = [loadNodeROMIds;modelDoFIDs];

Find the relevant components of rom.K and rom.M for time integration.

Kconstrained = rom.K(odeDoFs,odeDoFs);
Mconstrained = rom.M(odeDoFs,odeDoFs);
numODE = numel(odeDoFs);

Now you have a second-order system of ODEs. To use ode15s, you must convert this system into a system of first-order ODEs by applying linearization. This type of a first-order system is twice the size of the second-order system.

Mode = [eye(numODE,numODE),   zeros(numODE,numODE); ...
        zeros(numODE,numODE), Mconstrained];
Kode = [zeros(numODE,numODE), -eye(numODE,numODE); ...
        Kconstrained,         zeros(numODE,numODE)];
Fode = zeros(2*numODE,1);

The specified concentrated force load in the full system is along the z-direction, which is the third DoF in the ODE system. Accounting for the linearization, obtain the first-order system to get the loaded ODE DoF.

loadODEDoF = numODE + 3;

Specify the mass matrix and the Jacobian for the ODE solver.

odeoptions = odeset;
odeoptions = odeset(odeoptions,"Jacobian",-Kode);
odeoptions = odeset(odeoptions,"Mass",Mode);

Specify zero initial conditions.

u0 = zeros(2*numODE,1);

Solve the reduced system by using ode15s and the helper function.

function f = CMSODEf(t,u,Kode,Fode,centerVertex)
Fode(centerVertex) = 10*sin(6000*t);
f = -Kode*u +Fode;
end

tlist = 0:0.00005:3E-3;
sol = ode15s(@(t,y) CMSODEf(t,y,Kode,Fode,loadODEDoF), ...
                                 tlist,u0,odeoptions);

Compute the values of the ODE variable and the time derivatives.

[displ,vel] = deval(sol,tlist);

Reconstruct Solution for Full Model

Knowing the solution in terms of the interface DoFs and modal DoFs, you can reconstruct the solution for the full model. The reconstructSolution function requires the displacement, velocity, and acceleration at all DoFs in rom. Create the complete solution vector, including the zero values at the fixed DoFs.

u = zeros(size(rom.K,1),numel(tlist));
ut = zeros(size(rom.K,1),numel(tlist));
utt = zeros(size(rom.K,1),numel(tlist));
u(odeDoFs,:) = displ(1:numODE,:);
ut(odeDoFs,:) = vel(1:numODE,:);
utt(odeDoFs,:) = vel(numODE+1:2*numODE,:);

Create a transient results object using this solution.

RTrom = reconstructSolution(rom,u,ut,utt,tlist);

Compute the displacement in the interior at the center of the beam using the reconstructed solution.

coordCenter = [0;0;0];
iDispRTrom = interpolateDisplacement(RTrom,coordCenter);
figure
plot(tlist,iDispRTrom.uz)
title("Z-Displacement at Geometric Center")