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:

paper_calculate_period( t , y , whichTime , whichEndTime)
```function [ period , deflectionMean , deflectionStd ] = paper_calculate_period( t , y , whichTime  , whichEndTime)
% PAPER_CALCULATE_PERIOD  Calculates period of a signal
%
% [ PERIOD DEFLECTIONMEAN DEFLECTIONSTD ] = PAPER_CALCULATE_PERIOD( t , y , whichTime )
% calculates the PERIOD of the (t,y) signal based on the time at which the
% WHICHTIMEth peak occurs. The mean value and standard deviation of the
% steady-state deflection (DEFLECTIONMEAN and DEFLECTIONSTD, respectively) are
% calculated based on the four cycles that begin at the WHICHTIMEth peak.
%
% [ PERIOD DEFLECTIONMEAN DEFLECTIONSTD ] = PAPER_CALCULATE_PERIOD( ... , whichEndTime )
% calculates the deflection based on the first four cycles of the last
% WHICHTIME peaks right before WHICHENDTIME

% Copyright 2006, The MathWorks, Inc.

% Set the following flag to true to plot the signal and the limits used for the
% period and steady-state calculations
PLOT_ANALYSIS = false;

if PLOT_ANALYSIS
figure; plot( t , y , 'b-' );hold on;title('Period Analysis');
ybot=min(y);ytop=max(y);ydiff=(ytop-ybot)/10;
ybotpd=ybot-2*ydiff;ytoppd=ytop+2*ydiff;
ybot=ybot-ydiff;ytop=ytop+ydiff;
end

if nargin<4
whichEndTime = [];
end

dy = diff(y);
fdy = dy>0;
dfdy = diff(fdy);
wdfdy = find( dfdy==-1);

%
% Extract period
%
Twf=t( wdfdy( whichTime ) );
Twb=t( wdfdy( 1 ) );
period = (Twf - Twb) / ( whichTime - 1 );
if PLOT_ANALYSIS
plot( [Twf Twf] , [ybotpd ytoppd] , 'g-' );
plot( [Twb Twb] , [ybotpd ytoppd] , 'g-' );
end

%
% Extract steady-state deflection
%
if isempty( whichEndTime )
endIdx = length ( wdfdy );
else

% Find the peak that occurs on or before  whichEndTime
[ dummy , timeIdx ] = min( abs( t - whichEndTime ) );
whichEndPeaksIdx = find( wdfdy < timeIdx );
endIdx = whichEndPeaksIdx( end );

end

% Get peaks near the end
wdfdy( endIdx-whichTime+4 )
wdfdy( endIdx-whichTime+3 )
wdfdy( endIdx-whichTime+2 )
wdfdy( endIdx-whichTime+1 )
wdfdy( endIdx-whichTime   )
];

% Average from approximately the midpoint of the oscillation
steadyStateZeroIdx = round( diff( steadyStatePeakIdx )/4 ) + ...
for j=1:length( steadyStateZeroIdx )-1
deflection( j ) = mean( y( steadyStateZeroIdx( j+1 ) : ...
steadyStateZeroIdx( j )  ) ) ;
if PLOT_ANALYSIS
plot( [ t(steadyStateZeroIdx( j ))  t(steadyStateZeroIdx( j )) ] , ...
[ybot ytop] , 'r-' );
end
end

deflectionMean = mean( deflection );
deflectionStd  = std(  deflection );

if PLOT_ANALYSIS
plot( [ t(steadyStateZeroIdx( end ))  t(steadyStateZeroIdx( end )) ] , ...
[ybot ytop] , 'r:' );
keyboard;
end
```