3. Linear Transient Analysis of Beam

A simply supported beam is loaded by a central load P applied at time t=0 (see also input file beam.mdl). The model has been chosen such that the response under the given loads is linear, the displacements of the order 0.004 m being smaller than the thickness of the beam (0.02 m), and the beam thickness to beam length ratio is small, thus rendering the influence of the shear strain small. A theoretical transient solution for this configuration is available, formulated with the Euler beam theory[29]:

DY t = P 0 L 3 n 4 E I i = 1,3,5, ... 1- cos ω i t i 4

where i is the odd (symmetric around y=L/2) mode number, E the modulus of elasticity, I the moment of inertia around the bending axis, m the beam mass per length unit, and L the length of the beam, and F t = P 0

The force applied at t=0 and kept constant. The beam as displayed below is modeled with 40 B2.S.RS beam elements, see MDL input file beam.mdl.

Beam model.

Figure 74. Beam model.


The test contains the following files:

MDL file beam.mdl
Viewer scripts view.py plots a displacement-time graph.
anim.py plots the deformed shape as a function of time.
Additional hfunc.py contains the time history function for running the b2modal analysis (see Makefile).
analytical.py calculates an analytical solution (see Makefile).

A free vibration analysis gives an idea about the frequencies to be expected, at the same time checking the free vibration analysis solution and preparing the input for the modal analysis. Theoretical free vibration modes and their corresponding circular frequencies ω are available with

ω i = i 2 π 2 E I m L 4

For transient analysis, the force F is set to 1.0 at t=0 and kept constant. The corresponding case is defined by the natural boundary condition nbc 1

nbc
   set 1
      value 1000. dof 2 node 21
   end
end

and the case 1

cases
   case 1
      nbc 1
      ebc 1
   end
end

ebc 1 referring to the zero value constraints (see input file beam.mdl). The free vibration analysis of the model is defined by case 2 of the cases command:

cases
   case 2
      analysis free_vibration
      nmodes 10
      ebc 1
   end
end

The free vibration analysis solution is obtained by explicitly requiring case 2 in the adir command:

adir
   case 2
end

Eigenvalues are obtained by executing b2000++ and by printing then with the b2print_modes tool:

Mode Eigenvalue     Frequency      Omega          Stiffness      Mass          
   1 +4.2134718e+03 +1.0330948e+01 +6.4911261e+01 +5.4405341e+06 +1.2912236e+03
   2 +6.1166480e+04 +3.9361975e+01 +2.4731858e+02 +5.5698495e+07 +9.1060487e+02
   3 +2.6892823e+05 +8.2535033e+01 +5.1858291e+02 +1.2587087e+08 +4.6804631e+02
   4 +7.2025449e+05 +1.3507131e+02 +8.4867809e+02 +2.2402548e+08 +3.1103657e+02
   5 +8.2915332e+05 +1.4492308e+02 +9.1057856e+02 +1.0361753e+09 +1.2496788e+03

Note that the theoretical circular eigenfrequencies are derived from an Euler beam model. Thus, the shear, which has a certain influence for the present model, is not taken into account and explains small differences between theoretical and calculated circular eigenfrequencies.

With circular frequencies available from the free vibration analysis the B2000++ transient solver and the modal analysis tool b2modal is executed with Rayleigh damping values 0, 0.1, and 1.0 (critical damping). Loading the structure with a step function of 1.0 at t=0, the transient analysis eventually converges to the static solution, given a damping coefficient > 0. The modal analysis test selects the first 3 eigenvalues.

While the modal solver takes the damping coefficients of each mode ζ for the model from the argument list, the transient solver requires the Rayleigh damping coefficients α and β to be defined (see K. J. Bathe, Finite Element Procedures; Prentice-Hall):

α = 2 ω 1 ω 2 ( ω 1 ζ 2 ω 2 ζ 1 ) ω 1 2 ω 2 2

β = 2 ( ω 1 ζ 1 ω 2 ζ 2 ) ω 1 2 ω 2 2

ω and ζ are the circular frequencies and the damping percentage, respectively. A linear transient analysis is defined in case 3, which is split in the effective case definition (case 3) and the stage definition (stage 31). Adding a stage is necessary for specifying the time range (time_step):

   case 3
      title 'Transient analysis'
      stage 31 time_step (T)
   end

stage 31 integrates from time 0.0 to T, T being a parametrized value (see file beam.mdl).

The stage in its turn defines conditions of the stage:

   stage 31
      nbc 1 sfunction '1.'
      ebc 1
      multistep_integration_order 2
      step_size_init (DT)
      step_size_min (DT)
      step_size_max (DT)
   end

nbc 1 sfunction '1.' includes the nbc set 1, scaling it to 1.0 for the whole stage. ebc 1 includes the nbc set 1, keeping it constant during the stage (default). multistep_integration_order selects the time integration scheme, here of order 2. step_size_init defines the initial time step or increment, step_size_min the minimum time step or increment, and step_size_max the maximum time step or time increment.

Result response curves for the mid-span y displacement DY are displayed in the figures below:

Time-displacement response, no damping.

Figure 75. Time-displacement response, no damping.


Time-displacement response, no damping.

Figure 76. Time-displacement response, no damping.


Time-displacement response, no damping.

Figure 77. Time-displacement response, no damping.




[29] R.W. Clough, J. Penzien; Dynamics of Structures; McGraw-Hill (1975); ISBN 0-07-011392-0