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.m
```%% Script parameters: change as needed

% Copyright 2006, The MathWorks, Inc.

%
% Directory in which to find the data files
%

%
% 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 ) );

%
% Trim the mesh if requested
%
if doTrimMesh
trimMeshNodes = (theMesh( : , 2 ) == rightEdgeXCoord );
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) has %i points.\n' , rightEdgeXCoord , 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);
rightMeshIndeces = find( rightMeshNodes );
rightEdgeMeshNum = sum( rightMeshNodes );
rightEdgeMesh = theMesh( rightMeshIndeces , : );
subplot(3,1,3);plot3( rightEdgeMesh(:,2), rightEdgeMesh(:,3), rightEdgeMesh(:,4) , 'b.');
title('Right Edge Mesh');xlabel('x');ylabel('y');zlabel('z');

%
% Set up the input nodes and their scaling magnitudes
%
if distributeinput
rightEdgeDistributedForce = forcescale / rightEdgeMeshNum;
whichInput = rightMeshIndeces;
else
rightEdgeDistributedForce = forcescale;
whichInput = rightMeshIndeces( 1 );
end

systemInputs = zeros( [ 1 trimMeshNum ] );
systemInputs( whichInput ) = rightEdgeDistributedForce;

%
% Set up the output nodes
%
whichOutput = rightMeshIndeces( 1 );
outputCoordinate = rightEdgeMesh( whichOutput , : )
systemOutputs = zeros( [ 1 trimMeshNum ] );
systemOutputs( whichOutput ) = 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_pva , ss_cantilever_p ] = fea_to_statespace( ssData.modalMatrix , ...
ssData.eigenvalues , ssData.modalDamping , ...
ssData.inputNodes , ssData.outputNodes );

%
% Plot the step response
%
figure;
t=0:1e-3:5;
[ y_pva , t_pva , x_pva ] = step( ss_cantilever_pva , t );
paper_plot_position_data( t_pva , [ t_pva y_pva( : , 1 ) ] , 'State-space model: step response' , 5 , [] )
%plot( t_pva , y_pva( : , 1 ) );
set( gca , 'XLim' , [0 0.2] );
%xlabel( 'Time [sec]' , 'Interpreter' , 'latex');
%ylabel( 'Displacement [m]' , 'Interpreter' , 'latex');
%title( 'State-space model: step response' );
paperfig_dosave( 'paperfig_stepresponse' );

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