Code covered by the BSD License

# Modeling Flexible Bodies in SimMechanics

### Dallas Kennedy (view profile)

08 May 2006 (Updated )

Technical paper and examples on modeling flexibility in SimMechanics.

### Editor's Notes:

cosmos_cantilever_twoside.m
```%% Script parameters: change as needed

% Copyright 2006, The MathWorks, Inc.

%
% Directory in which to find the data files
%
if ispc
end

%
% Name prefix for the data files
%
namePrefix='cantilever_part';

%
% In this example, we use uniform damping (the same for all the modes).
% This is the coefficient that multiplies all the omegas in the damping
% matrix and is equal to 2*xi.
%
uniformDampingRatio = .15;

%
% Whether we should attempt to trim the mesh when obtaining the modes, so
% that the trimmed mesh only contains the input and output nodes. When the
% FEA mesh is very fine, this can make the MATLAB data structures smaller
% by throwing out data we don't use.
%
doTrimMesh = true;

%
% A scale factor for the force (e.g., to convert between systems of units).
% In this example, both the FEA model and the SimMechanics model use MKS
% units, so this is unity.
%
forcescale=1;

%
% Should we model the input as being distributed along the entire right
% edge of the beam or just at the first point we find on the right edge?
%
distributeinput=true;

%
% Which modes to include on the analysis
%
whichModes = [ 1 ];

%% Obtaining the FEA data

%
% Get the original mesh
%
[ md , theMesh ] = cosmos2m( datadir , namePrefix , true );
origMesh = theMesh;
figure(1);
subplot(3,1,1);plot3( theMesh(:,2), theMesh(:,3), theMesh(:,4) , 'b.');
title('Original Mesh'); xlabel('x');ylabel('y');zlabel('z');
fprintf( 1 , '==>Original mesh has %i points\n' , length( theMesh) );

% Select the nodes of particular interest
rightEdgeXCoord = max( theMesh( : , 2 ) );
leftEdgeXCoord = min( theMesh( : , 2 ) );

%
% Trim the mesh if requested
%
rightMeshNodes = ( theMesh( : , 2 ) == rightEdgeXCoord ) ;
leftMeshNodes  = ( theMesh( : , 2 ) == leftEdgeXCoord ) ;
if doTrimMesh
trimMeshNodes = (  rightMeshNodes | leftMeshNodes );
else
trimMeshNodes = ones( [ 1 length( theMesh ) ] );
end

trimMeshIndeces = find( trimMeshNodes );
trimMesh = theMesh( trimMeshIndeces , : );

trimMeshNum = length( trimMeshIndeces );

%
% Process the FEA data at the selected nodes only
%
[ md , theMesh ] = cosmos2m( datadir , namePrefix , [] , trimMeshIndeces , whichModes );
fprintf( 1 , '==>New mesh (x==%f,%f) has %i points.\n' , rightEdgeXCoord , leftEdgeXCoord , trimMeshNum );
subplot(3,1,2);plot3( theMesh(:,2), theMesh(:,3), theMesh(:,4) , 'b.');
title('Trimmed Mesh');xlabel('x');ylabel('y');zlabel('z');

%% Process the FEA data for the state-space model

%
% Make sure to identify the nodes at the right edge (necessary if we did
% not trim the mesh earlier )
%

rightMeshNodes = ( theMesh( : , 2 ) == rightEdgeXCoord ) ;
leftMeshNodes  = ( theMesh( : , 2 ) == leftEdgeXCoord ) ;
interestMeshNodes = (  rightMeshNodes | leftMeshNodes );

interestMeshIndeces = find( interestMeshNodes );
interestEdgeMeshNum = sum( interestMeshNodes );
interestEdgeMesh    = theMesh( interestMeshIndeces , : );

rightMeshIndeces = find( rightMeshNodes );
rightEdgeMeshNum = sum( rightMeshNodes );
rightEdgeMesh    = theMesh( rightMeshIndeces , : );

leftMeshIndeces = find( leftMeshNodes );
leftEdgeMeshNum = sum( leftMeshNodes );
leftEdgeMesh    = theMesh( leftMeshIndeces , : );

subplot(3,1,3);plot3( interestEdgeMesh(:,2), interestEdgeMesh(:,3), interestEdgeMesh(:,4) , 'b.');
title('Interest Edge Mesh');xlabel('x');ylabel('y');zlabel('z');

%
% Set up the input nodes and their scaling magnitudes
%
systemInputs = zeros( [ 2 trimMeshNum ] );
if distributeinput
rightEdgeDistributedForce = forcescale / rightEdgeMeshNum;
leftEdgeDistributedForce = forcescale / leftEdgeMeshNum;
systemInputs( 1 , rightMeshIndeces  ) = rightEdgeDistributedForce;
systemInputs( 2 , leftMeshIndeces   ) = leftEdgeDistributedForce;
else
rightEdgeDistributedForce = forcescale;
leftEdgeDistributedForce = forcescale;
systemInputs( 1 , rightMeshIndeces( 1 )  ) = rightEdgeDistributedForce;
systemInputs( 2 , leftMeshIndeces( 1 )   ) = leftEdgeDistributedForce;
end

%
% Set up the output nodes
%
whichOutput = [ rightMeshIndeces( 1 ) leftMeshIndeces( 1 ) ];
outputCoordinate = interestEdgeMesh( whichOutput , : );
systemOutputs = zeros( [ 2 trimMeshNum ] );
systemOutputs( 1 , whichOutput(1) ) = 1;
systemOutputs( 2 , whichOutput(2) ) = 1;

%
% Encapsulate the state-space data
%
ssData.modeList = 1:md.num_modes;
ssData.nodeList = 1:trimMeshNum;
ssData.modalMatrix = md.nodal_displacements_Y( ssData.nodeList , ssData.modeList );
ssData.eigenvalues = md.frequency_angular';
ssData.modalDamping = uniformDampingRatio * ones( size( ssData.modeList ) );
ssData.inputNodes = systemInputs';
ssData.outputNodes = systemOutputs;

%% Create the state-space system
[ ss_cantilever_twoside_pva , ss_cantilever_twoside_p ] = fea_to_statespace( ssData.modalMatrix , ...
ssData.eigenvalues , ssData.modalDamping , ...
ssData.inputNodes , ssData.outputNodes );

%
% Save the state-space system
%
save ss_cantilever_twoside_pva ss_cantilever_twoside_pva
disp( '==> Saved ss_cantilever_pva' )```