By Victor Chudnovsky, MathWorks, Dallas Kennedy, MathWorks, Arnav Mukherjee, MathWorks, and Jeff Wendlandt, MathWorks
SimMechanics and Simulink combine to form an efficient tool for simulating rigid-body mechanical systems, especially in control systems applications. These products also enable you to simulate flexible body motion—a frequent requirement in aerospace, automotive, civil engineering, and other industries. Because the SimMechanics Body block represents only a rigid body, you must first transform your continuous flexible body into an approximate form that you can represent in Simulink. SimMechanics can then simulate flexibility with a fidelity that is sufficient for many applications, such as control design.
In this article we use SimMechanics to apply the two most common flexible body approximation methods to modeling beams: the lumped-parameter approximation and the state space/frequency response method using finite element analysis (FEA) results. Both methods assume that beam deflection is small and in the linear regime.
The lumped parameter approach approximates a flexible body as a set of rigid bodies coupled with springs and dampers. It can be implemented by a chain of alternating bodies and joints. The springs and dampers act on the bodies or the joints. The spring stiffness and damping coefficients are functions of the material properties and the geometry of the flexible elements. While lumped parameter methods can be extended to more complicated systems, they are best suited for systems with chain-like geometries. It is easier to model bodies with multidimensional geometries using other approaches, such as FEA.
In this section we apply the lumped parameter method first to an abstract beam and then to a flexing cantilever. A beam is represented by a series of generalized beam elements (GBEs), each of which is a body-joint-body combination. This formulation actuates and senses the joints, taking advantage of the relative degrees of freedom (DoFs) inherent in SimMechanics.
The lumped parameter method involves five steps:
Simulink and SimMechanics simplify the task of deriving lumped parameter approximations by enabling you to:
To analyze the GBE motion and forces, define X and F to be generalized coordinates and forces at one end of the GBE. Assume that the other end is temporarily fixed. Let x represent the joint motion. Then
X = g(x) , dX = dg(x)·dx = [J]·dx , where [J] = dg/dx
(1)
The generalized force F at the tip and generalized force f at the joint can be expressed with a generalized stiffness matrix [K] or [k]:
F = [K]·dX , f = [k]·dx
(2)
The Jacobian [J] relates the two matrices, such that f = [J]^{T}·[K]·[J]·dx. Thus:
[k] = [J]^{T}·[K]·[J]
(2a)
Now we apply this technique to a cantilever or beam bending in a plane. In our example, one revolute represents the bending of a single GBE relative to its neighbor.
P and T are the force and moment (torque), respectively, at the end of the GBE. The moment τ at the joint can be expressed from equation (2a) as
τ = k dθ , where k = (E ·I _{Az})/L (3)
in terms of the area moment of inertia I_{Az} (second moment of area of the beam’s yz cross-section) and Young’s elastic modulus E. The rotational joint n executes damped oscillations according to the normalized moment equation
(4)
ω_{0}^{2} is the ratio of the effective spring constant k and the mass moment of inertia I_{Mz}, both assumed to be the same for all GBEs. The damping coefficient 2ξ_{0}ω_{0} that multiplies the velocity is a quasiempirical value that accounts for energy lost to viscoelastic effects.
The beam is made of an aluminum alloy of density 2700 kg/m^{3} and Young’s modulus 6.9x10^{10} N/m^{2}. We apply a load of 196.2 N (equivalent to the weight of 20 kg) to the tip.
Calculation yields these values, neglecting the beam weight.
Mass (analytic) | 1.35 kg |
Equilibrium deflection at tip (analytic) | -3.56 mm |
Equilibrium deflection at tip (lumped GBEs) | -3.54 mm |
Vibration frequency (lumped GBEs) | 63.7 Hz |
Now we apply the lumped parameter method with ten GBEs to the same cantilever, mounted to the wall with a torsional spring and including the beam weight. Again, we apply the 20 kg-mass-equivalent load to the tip. The beam bends by this load as well as by its own weight. (Bending under its own weight is a rigid body mode of the beam.) The addition of the beam weight reduces the vibration frequency considerably. This system and approximation use a model available from MATLAB Central.
Base spring stiffness K_{s} | 10,000 N/m |
Equilibrium deflection at tip (lumped GBEs) | -3.55 mm |
Vibration frequency (lumped GBEs) | 6.82 Hz |
The lumped parameter method outlined here uses purely local forces as functions of local deflections. This approximation is useful for estimating the gross behavior and endpoint deflection of the beam but not for modeling local dynamics along the beam or the beam’s vibrational modes. The restoring moment generated by bending the cantilever is proportional to its curvature. This method approximates the curvature by the local difference in slopes rather than by the correct nonlocal difference in slopes between neighboring GBEs, which cannot be represented in a series of independent elements. Thus, while the lumped parameter method is a convenient way to represent a flexible body, it does not, in general, yield correct results.
A more powerful approach to modeling flexible bodies in SimMechanics is to first use an FEA program to obtain the frequency or modal response of the deformed bodies and then to incorporate this response into a SimMechanics model by superimposing flexible body deflection on rigid body motion. This approach involves three steps:
The most direct way to implement the black box in Simulink is to use the State Space or LTI System block. Either block implements the flexible body dynamics in a state-space model and uses the results to actuate the motion of massless bodies to which other SimMechanics blocks can be connected. You can use this method to model multiple flexible bodies and complex systems comprising many flexible parts, if the deformations remain small enough for linear approximation..
The state-space model is based on the multi-DoF version of (4),
, with M, C, K, and F as mass, friction, spring, and force matrices acting on the joint DoFs q and system inputs u (for example, reaction forces). A proper choice of coordinate transformation Φ from the spatial DoFs q to the normal coordinates η = Φ·q yields the normalized form of the dynamical equations,
, and diagonalizes Ω, the stiffness matrix whose eigenvalues are the squared angular frequencies, ω_{n}^{2}. We obtain the η, Φ, and Ω values from the FEA program. Since there is no simple way to derive Γ from theory, we use approximations and empirical values instead. In the commonly used proportional damping ansatz, Γ is assumed to be diagonal, with each eigenvalue γ_{n} = ξ_{n}ω_{n}. A complete state space also includes system outputs y (for example, sensors), represented as linear functions of the DoFs and their velocities and accelerations:
We represent this collection of forced, damped, coupled linear harmonic oscillators in linear time-invariant (LTI) state-space form
with
where
and
Now we apply the FEA method to the same bending aluminum cantilever that we previously modeled using lumped parameters. Without spring-mounting at the wall or any beam weight, the lowest two vibration frequencies are 69.8 Hz and 434 Hz.
Consider this cantilever with the same spring loading and revolute joint at the wall and the same beam weight as before. Modeling the cantilever using only the lowest vibration mode in FEA yields a vibration frequency and equilibrium deflection of 8.76 Hz and -3.20 mm, respectively.
All FEA results were obtained from COSMOSWorks 2006 SP2.0, a commercial FEA software application, using default mesh parameters (global size 7.94 mm, tolerance 0.40 mm). The SimMechanics model is available from MATLAB Central.
With a reasonable number of discrete elements, the lumped parameter method gives accurate results for static bending, but the FEA method gives more accurate results for vibration modes and frequencies.
Body-loaded, spring-mounted cantilever | Vibration frequency (Hz) | Equilibrium deflection (mm) |
---|---|---|
Lumped parameters (10 GBEs) | 6.79 | -3.63 |
FEA/Frequency response | 8.76 | -3.20 |
Each method has advantages and disadvantages. The lumped parameter method is easy and quick to implement for flexible bodies with chain-like geometries but is more difficult and ambiguous, and hence less useful, for modeling bodies with higher-dimensional geometries. The simplification to independent GBEs also incorrectly represents the curvature bending moment. To correctly approximate the curvature moment you must apply the moment on each joint as a function of its own and neighboring joint deflections.
The FEA/state-space method is straightforward to implement for any body geometry and any number of DoFs. It is also more naturally suited to control design and analysis problems. Once you embed your FEA state space into a SimMechanics model, this method introduces algebraic loops as Simulink seeks a consistent solution for the acceleration of the Massless Bodies. (The acceleration depends on the reaction force, which is itself a function of the acceleration.) Algebraic loops can slow down the simulation and make it less accurate. You can break algebraic loops with Transfer Function blocks, with little or no effect on fidelity.
For a highly refined discretization and a large number of DoFs (many GBEs and joints in the lumped parameter case, many states in the FEA case) the dynamics become difficult to solve numerically.
Both methods can lead to dynamics controlled by a wide range of frequencies. A large ratio of the highest important frequency to the lowest makes the system technically "stiff" and difficult to simulate, meaning that the dynamics simultaneously reflect many different time scales. In the FEA method, you can select the modes and frequencies of interest, neglecting the others, to mitigate this difficulty.
A text-and-model archive, including a detailed technical paper and models referenced in this article, is available from MATLAB Central.
Published 2006