Spring-mass-damper system example

Octave script

ODE System definition

In this example a simple spring-mass-damper system is considered. This simple problem is used to validate different implementations of truss and frame elements submitted to dynamic loads.

spring-mass diagram

Analytic solution

The analytical solution is based on chapter 3 from Dynamics of Structures by Ray W. Clough and Joseph Penzien, Third Edition, 2003. The analytical solution of the problem is given by:

\[ u(t) = \left( A_c \cos( \omega_D t ) + B \sin( \omega_D t ) \right) e^{ -\xi \omega_N t } + G_1 \cos( \bar{\omega} t ) + G_2 \sin( \bar{\omega} t )\]

where the notation and the parameters interpretation can be seen in the cited reference.

A set of numerical parameters must be defined to compute the analytic solution.

We start as all models, clearing the workspace and adding the ONSAS path to the work path.

% clear workspace and add path
close all, clear all; addpath( genpath( [ pwd '/../../src'] ) );

The following numeric parameters are considered.

% scalar parameters for spring-mass system
k    = 39.47 ; % spring constant
c    = 2   ; % damping parameter
m    = 1     ; % mass of the system
p0   = 40    ; % amplitude of applied load
u0   = 0.1   ; % initial displacement
du0  = 0.0   ; % initial velocity

Then other parameters are computed:

omegaN  = sqrt( k / m )           ; % the natural frequency
xi      = c / m  / ( 2 * omegaN ) ;
freq    = omegaN / (2*pi)         ;
TN      = 2*pi / omegaN           ;
dtCrit  = TN / pi                 ;

The frequency of the sinusoidal external force is set as:

omegaBar = 4*omegaN ;

The analytic solution can be computed for specific cases as follows:

if (c == 0) && (p0 == 0) % free undamped solution
  myAnalyticFunc = @(t)   (   u0 * cos( omegaN * t )  ) ;
else                     % other cases solution
  beta   = omegaBar / omegaN ;  omegaD = omegaN * sqrt( 1-xi^2 ) ; %forced and damped
  G1 = (p0/k) * ( -2 * xi * beta   ) / ( ( 1 - beta^2 )^2 + ( 2 * xi * beta )^2 ) ;
  G2 = (p0/k) * (  1      - beta^2 ) / ( ( 1 - beta^2 )^2 + ( 2 * xi * beta )^2 ) ;
  
  Ac = u0 - G1 ;  B = (xi*omegaN*Ac - omegaBar*G2 ) / (omegaD) ;
  
  myAnalyticFunc = @(t) ...
     ( Ac * cos( omegaD   * t ) + B  * sin( omegaD   * t ) ) .* exp( -xi * omegaN * t ) ...
    + G1  * cos( omegaBar * t ) + G2 * sin( omegaBar * t ) ;
end

Numerical solutions

The analytic solution is used to validate two numerical solution approaches using different structural physical models, governed by the same ODE.

Numerical case 1: truss element model with Newmark method and lumped masses

In this case, a truss element is considered, as shown in the figure, with Young modulus, cross-section, area, mass, nodal damping and length corresponding to the parameters considered for the spring-mass-damper system

spring-mass diagram

The scalar parameters for the equivalent truss model are:

l   = 10                ;
A   = 0.2               ;
rho = m * 2 / ( A * l ) ;
E   = k * l /   A       ;

where the material of the truss was selected to set a mass $m$ at the node $2$.

assert( u0 < l, 'this analytical solution is not valid for this u0 and l0');

Materials

materials(1).hyperElasModel  = '1DrotEngStrain' ;
materials(1).hyperElasParams = [ E 0 ]          ;
materials(1).density         = rho              ;

Elements

In this case only 'node' and 'truss' elements are considered and the lumped inertial formulation is set for the truss element:

elements(1).elemType = 'node'                                 ;
elements(2).elemType = 'truss'                                ;
elements(2).elemCrossSecParams = {'circle', [sqrt(4*A/pi) ] } ;
elements(2).massMatType = 'lumped'                            ;

Boundary conditions

The node $1$ is fixed, so the boundary condition set is:

boundaryConds(1).imposDispDofs =  [ 1 3 5 ] ;
boundaryConds(1).imposDispVals =  [ 0 0 0 ] ;

The node $2$ allows the truss to move in $x$ so the boundary condition set is:

boundaryConds(2).imposDispDofs =  [ 3 5 ] ;
boundaryConds(2).imposDispVals =  [ 0 0 ] ;

ant the external load is added into the same boundary condition using:

boundaryConds(2).loadsCoordSys = 'global'                   ;
boundaryConds(2).loadsTimeFact = @(t) p0*sin( omegaBar*t )  ;
boundaryConds(2).loadsBaseVals = [ 1 0 0 0 0 0 ]            ;

Initial conditions

Initial displacement and velocity are set:

initialConds(1).nonHomogeneousUDofs    = 1   ;
initialConds(1).nonHomogeneousUVals    = u0  ;
initialConds(1).nonHomogeneousUdotDofs = 1   ;
initialConds(1).nonHomogeneousUdotVals = du0 ;

Analysis settings

The following parameters correspond to the iterative trapezoidal Newmark method with the following tolerances, time step, tolerances and final time

analysisSettings.methodName    = 'newmark' ;
analysisSettings.deltaT        =   0.005   ;
analysisSettings.finalTime     =   2.5*TN  ;
analysisSettings.stopTolDeltau =   1e-10    ;
analysisSettings.stopTolForces =   1e-10    ;
analysisSettings.stopTolIts    =   10      ;

OtherParams

The nodalDispDamping is added into the model using:

otherParams.nodalDispDamping =   c    ;

The name of the problem is:

otherParams.problemName = 'springMass_case1'     ;

mesh

Only two nodes are considered so the nodes matrix is:

mesh.nodesCoords = [  0  0  0 ; ...
                      l  0  0 ] ;
mesh.conecCell = { } ;
% The first node has no material, the first element of the _elements_ struct, which is `'node'` also the first boundary condition (fixed) and no initial condition is set.
mesh.conecCell{ 1, 1 } = [ 0 1 1 0   1   ] ;
% The second node has no material, the first element of the _elements_ struct, which is `'node'` also the second boundary condition (x disp free) and the first initial condition ($u_0$) is set.
mesh.conecCell{ 2, 1 } = [ 0 1 2 1   2   ] ;
% Only one element is considered with the first material and the second element setting
mesh.conecCell{ 3, 1 } = [ 1 2 0 0   1 2   ] ;

Execute ONSAS and save the results:

[matUsNewmark, loadFactorsMat] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;

Numerical case 2: truss model with nodal masses, using $\alpha$-HHT method and user loads function

Material

The nodalMass field allows to add lumped matrices to a node, since this field is used, then the equivalent $\rho$ of the material(1) aforementioned now is set to 0. Although an equal mass $m$ is considered for $u_x$ $u_y$ and $u_z$ at the node $2$, so:

materials(1).density   = 0       ;
materials(2).nodalMass = [m m m] ;

Boundary conditions

the boundary conditions struct is entirely re-written.

% repeat the BCs for node 1
boundaryConds = { } ;
boundaryConds(1).imposDispDofs =  [ 1 3 5 ] ;
boundaryConds(1).imposDispVals =  [ 0 0 0 ] ;
% repeat the BCs for node 2
boundaryConds(2).imposDispDofs =  [ 3 5 ] ;
boundaryConds(2).imposDispVals =  [ 0 0 ] ;

ant the external load is added into the same boundary condition using:

boundaryConds(3).userLoadsFilename = 'myLoadSpringMass' ;

where inside the function 'myLoadSpringMass' the external force vector of the structure with 12 = (2x6) entries is computed.

now the initial condition is added to the node $2$ with the second material:

mesh.conecCell{ 2, 1 } = [ 2 1 2 1   2  ] ;

The $\alpha_{HHT}$ method with $\alpha=0$ is equivalent to Newmark, this is employed to validate results of both methods, then:

analysisSettings.methodName    = 'alphaHHT' ;
analysisSettings.alphaHHT      =   0        ;
otherParams.problemName = 'springMass_case2'     ;

Execute ONSAS and save the results:

[matUsHHT, loadFactorsMat] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;

Numerical case 3: beam element model

d = l/10;
A   = pi * d^2 /  4 ;
Izz = pi * d^4 / 64 ; 
E   = k* l^3 / ( 3 * Izz ) ; % delta = P L3/(3EI)  =>  k = P/delta = 3EI/L3  => E = kL3/(3I)
rho = 2*m/(A*l) ; 
materials = {};
materials(1).hyperElasParams = [ E 0 ] ;
materials(1).density  = rho ;
materials(1).hyperElasModel  = 'linearElastic' ; % 1DrotEngStrain should work as well
elements = {} ;
elements(1).elemType = 'node' ;
elements(2).elemType = 'frame'; %and not truss
elements(2).massMatType = 'lumped' ;
elements(2).elemCrossSecParams{1,1} = 'circle' ;
elements(2).elemCrossSecParams{2,1} = [d] ;
elements(2).elemTypeParams = 0 ;
boundaryConds = {} ;
boundaryConds(1).imposDispDofs =  [ 1 2 3 4 5 6] ;
boundaryConds(1).imposDispVals =  [ 0 0 0 0 0 0] ;
boundaryConds(2).loadsCoordSys = 'global'                  ;
boundaryConds(2).loadsTimeFact = @(t) p0*sin( omegaBar*t )        ;
boundaryConds(2).loadsBaseVals = [0 0 1 0 0 0 ] ; %along Y axis

An initial displacements $u_0$ is set in $y$ direction:

initialConds = {} ;
initialConds(1).nonHomogeneousUDofs    = 3   ; 
initialConds(1).nonHomogeneousUVals    = u0  ;
initialConds(1).nonHomogeneousUdotDofs = 3   ;
initialConds(1).nonHomogeneousUdotVals = du0 ;
mesh.nodesCoords = [  0  0  0 ; ...
                      l  0  0 ] ;
mesh.conecCell = { } ;
mesh.conecCell{ 1, 1 } = [ 0 1 1 0   1   ] ;
mesh.conecCell{ 2, 1 } = [ 0 1 2 1   2   ] ;
mesh.conecCell{ 3, 1 } = [ 1 2 0 0   1 2   ] ;
otherParams.problemName = 'springMass_case3'     ;
analysisSettings.methodName    = 'newmark' ;
[matUsBending, loadFactorsMat] = ONSAS( materials, elements, boundaryConds, initialConds, mesh, analysisSettings, otherParams ) ;
valsBending = matUsBending(6+3,:) ;

Verification

The numerical displacements of the node $2$ is extracted for both study cases:

valsNewmark = matUsNewmark(6+1,:) ;
valsHHT     = matUsHHT(6+1,:)     ;

The analytical solution is evaluated:

times       = linspace( 0,analysisSettings.finalTime, size(matUsHHT,2) ) ;
valsAnaly   = myAnalyticFunc(times)                                                          ;

The boolean to validate the implementation is evaluated such as:

analyticCheckTolerance = 5e-2 ;
verifBooleanNewmark =  ( ( norm( valsAnaly  - valsNewmark ) / norm( valsAnaly ) ) <  analyticCheckTolerance ) ;
verifBooleanHHT     =  ( ( norm( valsAnaly  - valsHHT     ) / norm( valsAnaly ) ) <  analyticCheckTolerance ) ;
verifBooleanBending =  ( ( norm( valsBending- valsHHT     ) / norm( valsAnaly ) ) <  analyticCheckTolerance ) ;
verifBoolean        = verifBooleanHHT && verifBooleanNewmark && verifBooleanBending                                   ;

Plot verification

The control displacement $u(t)$ is plotted:

figure
hold on, grid on, spanPlot = 8 ; lw = 2.0 ; ms = 11 ; plotfontsize = 20 ;
plot(times, valsAnaly   ,'b-', 'linewidth', lw,'markersize',ms )
plot(times(1:spanPlot:end), valsNewmark(1:spanPlot:end) ,'ro', 'linewidth', lw,'markersize', ms )
plot(times(1:spanPlot:end), valsHHT(1:spanPlot:end)     ,'gs', 'linewidth', lw,'markersize', ms )
plot(times(1:spanPlot:end), valsBending(1:spanPlot:end)     ,'yx', 'linewidth', lw,'markersize', ms )
labx = xlabel('t [s]');   laby = ylabel('u(t) [m]') ;
legend( 'analytic', 'truss-Newmark', 'nodalMass-HHT', 'Beam-model', 'location','southeast')
title( sprintf('dt = %.3d, m = %d, c = %d, k = %d, p0 = %d', analysisSettings.deltaT, m, c, k, p0 ) )

set(gca, 'linewidth', 1.0, 'fontsize', plotfontsize )
set(labx, 'FontSize', plotfontsize); set(laby, 'FontSize', plotfontsize) ;
if exist('../../docs/src/assets/')==7
  % printing plot also to docs directory
  disp('printing plot also to docs directory')
  print('../../docs/src/assets/springMassCheckU.png','-dpng')
else
  print('output/springMassCheckU.png','-dpng')
end
plot check