In this tutorial, the capabilities of the VIBES Toolbox for Component Mode Synthesis are demonstrated. Component Mode Synthesis (CMS) is a Dynamic Substructuring method for modally reduced numerical models. This tutorial shows this method by assembling a building from repeated substructures and assessing its dynamic behavior.

## Load FE models

Start by loading the FE models of the substructures as `vibes.MCKModel`

objects. Each model comprises a mass matrix, a stiffness matrix, the nodes of the FE elements, and a DoF-mapping table. Two substructures are defined:

- substructure A: a side part of the building
- substructure B: a mid part of the building

Both substructures comprise two stories and can be combined to create a taller and wider building. The models also come with a set of FE element definitions for the floors and walls, which are used to show the mode shapes as faces of the `vibes.ui.Mesh`

objects.

% Load the substructure models mckA = vibes.load('VIBES,Examples,Datasets','mat','Building','substructureA.MCKModel.mat'); mckB = vibes.load('VIBES,Examples,Datasets','mat','Building','substructureB.MCKModel.mat');

Plot a mesh of both objects in the 3D viewer to see their shape and orientation. Note that faces are obtained from the field `Parameters.Faces`

in the objects.

V = vibes.figure3d(); V.clf(); mckA.plotMesh('Color','b'); mckB.plotMesh('Color','r'); % Set camera V.view(0);

## Rotating the models

Looking at the meshes, it appears that the models are oriented incorrectly. We can rotate both models around the corner of the building to place them horizontally using the `vibes.MCKModel.rotate`

command.

% Rotate about the X-axis (roll) to place the buildings on the ground (XY-plane) mckA.rotate([90 0 0]); mckB.rotate([90 0 0]); % Move the models in Y-direction mckA.move([0 6 0]); mckB.move([0 6 0]);

Plot the new models in the 3D viewer.

% Clear the 3D viewer V.clf(); % Plot the rotated models mckA.plotMesh('Color','b'); mckB.plotMesh('Color','r'); % Set camera V.view(0); V.view('persp'); V.view(-20,15); V.zoom(2); % Save camera cam1 = V.Camera;

## Duplicating the models

Calling the `vibes.MCKModel.move`

or `vibes.MCKModel.rotat`

e methods with an output variable will create a shifted or rotated duplicate of that model, leaving the original unchanged. We will use this syntax to create shifted copies of the two substructures and build the larger building.

% Add a side element to close the right side of the building mckC = mckA.move([25 0 0]); mckC.Name = 'Substructure C'; % Add another floor mckD = mckA.move([0 0 6]); mckD.Name = 'Substructure D'; mckE = mckB.move([0 0 6]); mckE.Name = 'Substructure E'; mckF = mckA.move([25 0 6]); mckF.Name = 'Substructure F';

Plot the additional models too to preview the whole assembly.

% Plot the complete building mckC.plotMesh('Color','b'); mckD.plotMesh('Color','b'); mckE.plotMesh('Color','r'); mckF.plotMesh('Color','b'); % Set camera V.view(0); V.view('persp'); V.view(-20,15); V.zoom(2); % Save camera cam2 = V.Camera;

## Hurty/Craig-Bampton reduction

To construct the complete building efficiently, we will perform a Hurty/Craig-Bampton model reduction on the two substructure models that make up the assembly. The H/CB method reduces the full-sized system matrices to reduced system matrices comprising boundary DoFs and fixed-interface vibration modes.

The coupling will take place at the boundary DoFs. To find these boundary DoFs, the complete set of DoFs is first passed through the `vibes.DoF.uniqueMatch`

method, which applies the default DoF-matching logic to get the “primal” set of DoFs to retain after coupling. The DoFs that are not unique can be found by taking the set difference of all DoFs and the unique ones, yielding the boundary DoFs. Finally, by taking the “DoF-matched” intersection of all DoFs with these boundary DoFs, the full set of boundary DoFs is found.

% Get all substructure DoFs uA = mckA.DoFs; uB = mckB.DoFs; uC = mckC.DoFs; uD = mckD.DoFs; uE = mckE.DoFs; uF = mckF.DoFs; % Stack to get the complete set of DoFs u = [uA; uB; uC; uD; uE; uF];

Use `vibes.DoF.uniqueMatch`

to find the DoFs that are unique. This applies DoF-matching on `vibes.DoF`

objects, so may take a while.

% Find the unique set of DoFs q = uniqueMatch(u); % Get the (redundant) coupling DoFs by taking the set difference of u and q u_2r = setdiff(u,q); % Get the complete set of boundary DoFs for all substructures u_2 = intersectMatch(u,u_2r,'SingleMatch',false);

As the substructures were copied and then moved, each substructure has the same node numbers as its original. We can use this to our advantage by reducing only substructure A at those nodes that are on the boundaries of A, C, D and F, and do the same for substructure B at the boundaries of B and E.

% Get the node numbers of the boundary DoFs nnr_2 = unique(u_2.getNodeNumbers()); % Get the boundary DoFs of A and B uA_2 = uA.select('NodeNumber',nnr_2); uB_2 = uB.select('NodeNumber',nnr_2); % Get boundary nodes from the models nodesA_2 = mckA.selectNodes('NodeNumber',nnr_2); nodesB_2 = mckB.selectNodes('NodeNumber',nnr_2); % Plot the boundary nodes nodesA_2.plot(vibes.plotStyle('ob6')); nodesB_2.plot(vibes.plotStyle('sr6')); % Make the meshes transparent alpha([V.Meshes(1:2).MeshPatch],0.5); alpha([V.Meshes(3:6).MeshPatch],0.1);

Reduce the full `vibes.MCKModel`

objects of just A and B to `vibes.CMSModel`

objects using the found boundary DoFs. We will be using a Hurty/Craig-Bampton reduction with 20 internal vibration modes for each substructure. Setting the `KeepTransformation`

property to true allows to visualize later on the fixed-interface vibration modes of the Hurty/Craig-Bampton model. However, it is computationally less efficient to do so, as these matrices may become considerably large.

% Define a number of fixed interface vibration modes n_r = 20; % Reduce the MCK models cmsA = mckA.reduce('HCB',uA_2,'nComponentModes',n_r,'KeepTransformation',true); cmsB = mckB.reduce('HCB',uB_2,'nComponentModes',n_r,'KeepTransformation',true);

Use the `vibes.CMSModel.spy`

method to show the structure of the reduced matrices. Notice the diagonal structure of the section of the matrices containing the fixed interface modes.

cmsA.spy();

Let us now inspect the modal basis of the reduced CMS models. Plotting a mesh from a `vibes.CMSModel`

object that has its reduction matrix stored in the `T`

field will automatically pass that information onto the mesh object. Press SPACE and CTRL-UP/DOWN to animate the fixed interface mode shapes.

% Clear 3D viewer V.clf(); % You can define mesh colormaps to show color gradients in the faces meshA = cmsA.plotMesh(); meshA.Animation.ColorMap = winter(); meshB = cmsB.plotMesh(); meshB.Animation.ColorMap = autumn(); % Activate the colormap for the faces meshA.Faces.setStyle('Color','fcn'); meshB.Faces.setStyle('Color','fcn'); % Set camera V.Camera = cam1; % Animate mode V.startAnimation(); pause(3.5); V.pauseAnimation();

## Duplicating the reduced models

To create the full building from the reduced models, we will duplicate and shift the newly created `vibes.CMSModel`

objects to their respective positions, similar as before.

cmsC = cmsA.move([25 0 0]); cmsC.Name = 'Substructure C'; cmsD = cmsA.move([0 0 6]); cmsD.Name = 'Substructure D'; cmsE = cmsB.move([0 0 6]); cmsE.Name = 'Substructure E'; cmsF = cmsA.move([25 0 6]); cmsF.Name = 'Substructure F';

We can assemble the six reduced substructure models together on their boundary nodes. To prevent double node numbers from showing up when coupling the substructures, we will reset the node numbers for each substructure before coupling. We will set them with offsets of 2000 per substructure, and also reset the names of the nodes to represent the node number.

cmsA.resetNodeNumbers(); cmsA.Nodes.resetNames(); cmsB.resetNodeNumbers(2000); cmsB.Nodes.resetNames(); cmsC.resetNodeNumbers(4000); cmsC.Nodes.resetNames(); cmsD.resetNodeNumbers(6000); cmsD.Nodes.resetNames(); cmsE.resetNodeNumbers(8000); cmsE.Nodes.resetNames(); cmsF.resetNodeNumbers(10000); cmsF.Nodes.resetNames();

Also plot the other parts to show the whole assembly.

cmsC.plotMesh(); cmsD.plotMesh(); cmsE.plotMesh(); cmsF.plotMesh(); % Also plot the boundary nodes nodes_2 = u_2r.getUniqueNodes(); nodes_2.plot('og3'); % Set camera V.Camera = cam2;

## Assembling the reduced models

The assembly of the substructures is performed by the `vibes.mixin.PhysicalDomain.couple`

method, which automatically matches the set of corresponding DoFs. Note that this only assembles the system matrices of the substructures, but does not yet give access to the coupled-system response.

[cms,coupleOpts] = couple(cmsA,cmsB,cmsC,cmsD,cmsE,cmsF); cms.Name = 'Assembled system';

The “couple” method returns a `vibes.operation.NumericalCouplingOperation`

object that allows you to inspect the coupling procedure. For instance, you can visualize the result of the automatic DoF-matching procedure in a matrix view.

vibes.figure('Spy'); clf(); coupleOpts.DoFMatching.spy();

The faces of the coupled system are simply the faces of all substructures combined.

cms.Parameters.Faces = [cmsA.Parameters.Faces; cmsB.Parameters.Faces;.... cmsC.Parameters.Faces; cmsD.Parameters.Faces;... cmsE.Parameters.Faces; cmsF.Parameters.Faces];

Using the spy method on the coupled system shows the overlapping of the matrices at the boundary DoFs, while the fixed interface modes are simply concatenated.

cms.spy();

Clear 3D viewer.

V.clf(); % Plot a mesh of the assembled system mesh = cms.plotMesh(); % Set camera V.Camera = cam2;

Since the mesh is quite large, automatically scaling the mode shapes will give a distorted image of the actual displacement of the nodes. Instead, we will turn off auto-scaling and scale up every mode by the same amount.

mesh.Animation.AutoScale = false; mesh.Animation.Scale = 10; % Animate mode mesh.start(); pause(3.5); mesh.pause();

## Solving the assembled vibration modes

The resulting CMS model after assembly consists of just under 3500 DoFs. A modal solution can be computed from the CMS model using` vibes.CMSModel.toModalModel`

.

We would like to constrain the building to the ground. To achieve this, we will select all boundary DoFs in the coupled system with a Z-coordinate below 0.1m. These DoFs will be constrained in the computation of the modal solution. Constraints can be provided using the optional parameter `ConstrainedDoFs`

.

% First, select the boundary DoFs to constrain u_constr = cms.BoundaryDoFs.select('Z','<',0.1); % Compute the constrained modal solution for the first 50 modes mm = cms.toModalModel(50,'Symmetrize','avg','ConstrainedDoFs',u_constr);

The modal solution we just calculated consists of the same DoF set as the CMS model it was calculated from. To obtain a solution that describes the displacements of the complete node-set, we like to expand the solution to the original DoFs. For this, we can use the reduction matrix, which is stored in the T field of the `vibes.CMSModel`

class.

% Select the appropriate transformation matrix from the 'T' field T = cms.T.select('Name','Reduction'); % Multiply the modal model with the transformation matrix mm_full = T * mm;

Plotting a mesh of this full modal model allows you to visualize mode shapes of the complete, coupled system.

% Clear 3D viewer V.clf(); % Plot a mesh of the assembled system mesh = mm_full.plotMesh(); % Plot the constrained nodes u_constr.plotNodes('xy4'); % Set camera V.Camera = cam2; % Animate mode mesh.start(); pause(3.5); mesh.pause();

## Calculating FRFs

Now that we have a complete model of the coupled substructures, it is very easy to calculate the Frequency Response at the boundary DoFs of the system. This can be done by matrix-inversion of the reduced dynamic stiffness.

% First, add some damping cms.addDamping('Rayleigh',1e-1,1e-6);

We will select three of the boundary DoFs as inputs and outputs for our FRF matrix: a Z-direction on floor 1 to an X-direction on a wall on the same side of the building and a Z-direction on the other side of the building.

% Get boundary nodes bnodes = cms.BoundaryNodes(); n1 = bnodes.selectNearest([0 2 3]); n2 = bnodes.selectNearest([5 2 6]); n3 = bnodes.selectNearest([35 2 9]); dofs(1) = cms.BoundaryDoFs.select('Node',n1,'Direction','==',[1 0 0],'first'); dofs(2) = cms.BoundaryDoFs.select('Node',n2,'Direction','==',[0 0 1],'first'); dofs(3) = cms.BoundaryDoFs.select('Node',n3,'Direction','==',[0 0 1],'first'); % Calculate FRFs for displacement over force between 1 and 15 Hz. frftype = 0; freq = 1:15; Y = cms.toFRFMatrix(dofs,dofs,freq,frftype,'ConstrainedDoFs',u_constr);

Plot the generated FRF.

% Open a new figure vibes.figure('FRF'); clf % Create a plot Y.plot([2 3],1,'label','<ch.label>');

Indicate the channels in the 3D viewer.

V.stopAnimation(); Y.Channels([2 3]).plot('Length',2); Y.RefChannels(1).plot('Length',2); % Make building transparent alpha(0.1)