Longitudinal bars, torsional rods and strings are governed by the wave equation, which can be written under the following form: \[ m \ddot y(x, t) - D \frac{\partial^2 y(x,t)}{\partial x^2} = p(x, t), \] where:
\(y(x, t)\): Kinematic data at location \(x\) and time \(t\)
\(p(x, t)\): External excitation term
\(m\): Linear inertia of the type
\(D\): Stiffness of the type
For a longitudinal bar:
\(y(x, t) = u(x, t)\)
\(p(x, t)\): Distributed longitudinal force [N/m]
\(m = \rho S\): Linear mass density [kg/m2]
\(D = E S\): Longitudinal stiffness [N]
For a torsional rod:
\(y(x, t) = \theta(x, t)\)
\(p(x, t)\): Distributed moment [N.m/m]
\(m = \rho I_G\): Linear rotational inertia [kg.m4]
\(D = G J_T\): Rotational stiffness [N.m2]
For a string:
\(y(x, y)\): Transverse displacement [m]
\(m\): Linear mass density [kg/m]
\(D\): Tension force [N]
Euler-Bernoulli beams are governed by the following equation of motion: \[ m\ddot v(x, t) + D\frac{\partial^4 v(x, t)}{\partial x^4} = p(x, t), \] where:
\(v(x, t)\): Transverse displacement [m]
\(p(x, t)\): External excitation term [N/m]
\(m = \rho S\): Linear mass density [kg/m]
\(D = E I_z\): Bending stiffness [N.m2]
Data types
All the following data types are a subtype of the super type OneDtype.
Bar
Bar(L, S, E, ρ)Structure containing the data of a homogeneous and isotropic longitudinal bar
Constructor parameters
L: Length [m]
S: Cross-section area [m²]
E: Young's modulus [Pa]
ρ: Mass density [kg/m³]
Fields
L: Length [m]
m: Line mass [kg/m]
D: Stiffness coefficient [Pa]
Rod
Rod(L, I, J, G, ρ)Structure containing the data of a homogeneous and isotropic torsional bar
Constructor parameters
L: Length [m]
I: Second-moment of area [m⁴]
J: Torsion constant [m⁴]
G: Shear modulus [Pa]
ρ: Mass density [kg/m³]
Fields
L: Length [m]
m: Line mass [kg/m]
D: Stiffness coefficient [Pa]
Strings
Strings(L, S, D, ρ)Structure containing the data of a homogeneous and isotropic string
Constructor parameters
L: Length [m]
S: Cross-section area [m²]
D: Tension [N]
ρ: Mass density [kg/m³]
Fields
L: Length [m]
m: Linear mass density [kg/m]
D: Tension [N]
Beam
Beam(L, S, I, E, ρ)Structure containing the data of a homogeneous and isotropic bending beam
Constructor parameters
L: Length [m]
S: Cross-section area [m²]
I: Second moment of area [m⁴]
E: Young's modulus [Pa]
ρ: Density [kg/m³]
Fields
L: Length [m]
M: Linear mass density [kg/m]
D: Bending stiffness [N.m²]
Related functions
modefreq
modefreq(model::Bar, fmax, bc = :CC)
modefreq(model::Rod, fmax, bc = :CC)
modefreq(model::Strings, fmax, bc = :CC)
modefreq(model::Beam, fmax, bc = :SS)Computes the natural frequencies of a longitudinal or torsional bar up to fmax
Inputs
model: Structure containing the bar data
fmax: Maximum frequency for calculating the mode shapes [Hz]
bc: Boundary conditions
For all OneDStructure
:CC: Clamped - Clamped
:CF: Clamped - Free
:FF: Free - Free
For beams
:SS: Simply Supported - Simply Supported
:SC: Simply Supported - Clamped
:SF: Simply Supported - Free
Outputs
ωn: Natural frequencies calculated up to ωmax = 2π*fmax [Hz]
kn: Vector of modal wavenumbers
modeshape
modeshape(model::Bar, kn, x, bc = :CC)
modeshape(model::Rod, kn, x, bc = :CC)
modeshape(model::Strings, kn, x, bc = :CC)
modeshape(model::Beam, kn, x, bc = :SS)Computes the mass-normalized mode shapes of a longitudinal or torsional bar
Inputs
model: Structure containing the bar data
kn: Array of modal wavenumbers
x: Coordinates of calculation points of the mode shapes
bc: Boundary conditions
For all OneDStructure
:CC: Clamped - Clamped
:CF: Clamped - Free
:FF: Free - Free
For beams
:SS: Simply Supported - Simply Supported
:SC: Simply Supported - Clamped
:SF: Simply Supported - Free
Output
ϕ: Mass-normalized mode shapes
# Dimensions
L = 1.
d = 3e-2
# Section features
S = π*d^2/4
Iz = π*d^4/64
IG = 2Iz
J = IG
# Tension for string
T = 100.
# Material
E = 2.1e11
ν = 0.33
G = E/(1 - 2*ν)
ρ = 7800.
# Computation parameters
fmax = 2000.
x = [0.1, 0.9]
# Initialization of the data types
bar = Bar(L, S, E, ρ)
rod = Rod(L, IG, J, G, ρ)
strings = Strings(L, S, T, ρ)
beam = Beam(L, S, Iz, E, ρ)
# Computation of the natural frequencies
ωn, kn = modefreq(bar, fmax)
# Computation of the corresponding mode shapes
ϕn = modeshape(bar, kn, x, :CC)Rectangular membranes are governed by the following equation of motion: \[ m \ddot w(x, y ,t) + D\Delta w(x, y, t) = p(x, y, t), \] where:
\(w(x, y, t)\): Transverse displacement [m] at point \((x, y)\) and time \(t\)
\(m\): Surface mass [kg/m]
\(D = \tau\): Linear tension [N/m]
\(\Delta = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\): Laplacian operator
Rectangular Kirchhoff-Love plates are governed by the following equation of motion: \[ m \ddot w(x, y, t) + D \Delta^2 w(x, y, t) = p(x, y, t), \] where:
\(w(x, y, t)\): Transverse displacement [m] at point \((x, y)\) and time \(t\)
\(m = \rho h\): Surface mass [kg/m2]
\(D = \frac{Eh^3}{12(1 - \nu^2)}\): Bending stiffness [N.m]
\(\Delta^2 = \frac{\partial^4}{\partial x^4} + 2\frac{\partial^4}{\partial x^2 \partial y^2} + \frac{\partial^4}{\partial y^4}\): Bilaplacian operator
Data type
All the following data types are a subtype of the super type TwoDStructure.
Membrane
Membrane(L, b, m, D)Structure containing the data of a homogeneous and isotropic rectangular membrane
Fields
L: Length [m]
b: Width [m]
m: Surface mass [kg/m²]
D: Tension per unit length [N/m]
Plate
Plate(L, b, h, E, ρ, ν)Structure containing the data of a homogeneous and isotropic bending plate
Constructor parameters
L: Length [m]
b: Width [m]
E: Young's modulus [Pa]
ρ: Density [kg/m³]
ν: Poisson's ratio
Fields
L: Length [m]
b: Width [m]
m: Surface mass [kg/m²]
D: Bending stiffness [N.m]
Related functions
modefreq
modefreq(model::Plate, fmax)
modefreq(model::Membrane, fmax)Computes the natural frequencies of a simply supported rectangular plate or a clamped rectangular membrane up to fmax
Inputs
model: Structure containing the data related to the plate
fmax: Maximum frequency for calculating the modal shapes [Hz]
Outputs
ωmn: Natural frequencies calculated up to ωmax = 2π*fmax [Hz]
kmn: Matrix of modal wave numbers
modeshape
modeshape(model::Plate, kpq, x, y)
modeshape(model::Membrane, kpq, x, y)Computes the mass-normalized mode shapes of a simply supported rectangular plate or a clamped rectangular membrane
Inputs
model: Structure containing the data related to the structure
kmn: Matrix of modal wave numbers
(x, y): Coordinates of the points where the mode shapes are calculated
Output
ϕ: Mass-normalized mode shapes
# Dimensions
Lp = 0.6
bp = 0.4
hp = 1e-3
# Material parameters
E = 2.1e11
ρ = 7800.
ν = 0.33
# Computation parameters
fmax = 1000.
xp = [0.1, 0.5]
yp = [0.1, 0.3]
# Initialization of the data types
plate = Plate(Lp, bp, hp, E, ρ, ν)
# Computation of the natural frequencies
ωn, kn = modefreq(plate, fmax)
# Computation of the corresponding mode shapes
ϕn = modeshape(plate, kn, xp, yp)Single degree of freedom (Sdof) systems are classically composed of a mass \(m\), a stiffness \(k\) and a viscous damper \(c\) (see Figure 1).
Mathematically, their dynamic behavior is governed by the following normalized equation of motion : \[ \ddot x(t) + 2\xi\,\omega_0\, \dot x(t) + \omega_0^2 x(t) = \frac{F(t)}{m}. \] where \(F(t)\) can be either a base or an external excitation applied to the system.
The Sdof system can thus be defined by:
Sdof
Sdof(m, ω0, ξ)Structure containing the data of a sdof system
Constructor
m: Mass [kg]
f0: Natural frequency [Hz]
ξ: Damping ratio
Fields
m: Mass [kg]
ω0: Natural frequency [rad/s]
ξ: Damping ratio
# Definition of the structural parameters
m = 1.
f₀ = 10.
ξ = 0.01
# Initialization of Sdof
sdof = Sdof(m, f₀, ξ)StructuralVibration.jl considers Multi-degrees of freedom (Mdof) systems, which topology is presented in Figure 2. This choice has been made, because it allows modeling a large variety of possible configurations.
The dynamic behavior of such a system is governed by the following matrix system: \[ \mathbf{M} \ddot{\mathbf{x}}(t) + \mathbf{K}\mathbf{x}(t) = \mathbf{f}(t), \] where:
\(\mathbf{M} = \text{diag}(m_1, \dots, m_j, \dots, m_N)\) is the mass matrix.
\(\mathbf{K}\) is the stiffness matrix such that: \[ \mathbf{K} = \begin{bmatrix} k_1 & -k_1 & 0 & \ldots & 0 & 0 \\ -k_1 & k_1 + k_2 & -k_2 & \ddots & \vdots & \vdots \\ 0 & -k_2 & \ddots & \ddots & 0 & \vdots \\ \vdots & 0 & \ddots & \ddots & -k_{N-1} & 0 \\ \vdots & \vdots & \ddots & -k_{N-1} & k_{N-1} + k_N & -k_N \\ 0 & 0 & \ldots & 0 & -k_N & k_N \end{bmatrix}. \]
\(\mathbf{x}(t) = \left[x_1(t), \dots, x_j(t), \dots, x_N(t)\right]^\mathsf{T}\) is the displacement vector.
\(\mathbf{f}(t) = \left[F_1(t), \dots, F_j(t), \dots, F_N(t)\right]^\mathsf{T}\) is the external force vector.
Data types
Mdof
Mdof(k, m, c = Float64[])Structure containing the data for building a mdof system
Fields
k: Stiffness coefficients of the spring elements
m: Masses of the mdof system
c: Damping coefficients of the viscous dampers
If viscous dampers are defined, the damping matrix \(\mathbf{C}\) is consistent with the stiffness matrix \(\mathbf{K}\), meaning that: \[ \mathbf{C} = \begin{bmatrix} c_1 & -c_1 & 0 & \ldots & 0 & 0 \\ -c_1 & c_1 + c_2 & -c_2 & \ddots & \vdots & \vdots \\ 0 & -k_2 & \ddots & \ddots & 0 & \vdots \\ \vdots & 0 & \ddots & \ddots & -c_{N-1} & 0 \\ \vdots & \vdots & \ddots & -c_{N-1} & c_{N-1} + c_N & -c_N \\ 0 & 0 & \ldots & 0 & -c_N & c_N \end{bmatrix}. \]
MdofMesh
MdofMesh(Elt, constrained_dofs, free_dofs)Structure containing the data for building a mdof mesh
Constructor
model: Mdof model
bc: Boundary conditions
:CC: Clamped - Clamped
:CF: Clamped - Free
:FF: Free - Free
Fields
Elt: Element connectivity matrix
constrained_dofs: Constrained degrees of freedom
free_dofs: Free degrees of freedom
Related functions
assembly
assembly(model::Mdof)Assembly of the mass, stiffness and damping matrices of a mdof system
Input
model: Mdof model
Outputs
K: Stiffness matrix
M: Mass matrix
C: Damping matrix (if viscous dampers are defined in model)
apply_bc
apply_bc(A, mesh)Apply boundary conditions to a given matrix
Inputs
A: Matrix to apply the boundary conditions
mesh: Mesh of the system
Output
A_bc: Matrix with boundary conditions applied
eigenmode
eigenmode(K, M, n = size(K, 1))Computes the eigenmodes of a system defined by its mass and stiffness matrices.
Inputs
K: Stiffness matrix
M: Mass matrix
n: Number of modes to be keep in the modal basis
Outputs
ω: Vector of natural frequencies
Φ: Mass-normalized mode shapes
modal_matrices
modal_matrices(ωn, ξn)Computes the modal mass, stiffness, and damping matrices
Inputs
ωn: Vector of natural frequencies
ξn: Modal damping
Outputs
Kn: Generalized stiffness matrix
Mn: Generalized mass matrix (identity matrix, due to mass normalization)
Cn: Generalized damping matrix
modal_effective_mass
modal_effective_mass(M, ϕn, r)Computes the effective mass of a mode
Inputs
M: Mass matrix
ϕn: Mode shape
r: Influence vector (rigid body mode)
Outputs
meff: Modal effective mass
Note: The modeshapes are supposed to be mass-normalized
# Definition of the structural parameters
k_mdof = [1., 1.]
m_mdof = ones(3)
c_mdof = [0.1, 0.1]
# Initialization of Mdof
mdof = Mdof(k_mdof, m_mdof, c_mdof)
# Definition of a MdofMesh
mdof_mesh = MdofMesh(mdof, bc = :CF)
# System assembly
K_mdof, M_mdof, C_mdof = assembly(mdof)
# Apply boundary conditions (if any)
K_bc = apply_bc(K_mdof, mdof_mesh)
M_bc = apply_bc(M_mdof, mdof_mesh)
C_bc = apply_bc(C_mdof, mdof_mesh)
# Compute the eigenmodes of the systems
ωn, Φn = eigenmode(K_bc, M_bc)
# Computation of the modal matrices
Kmodal, Mmodal, Cmodal = modal_matrices(ωn, 0.01)
# Computation of modal effective mass
meff = modal_effective_mass(M_bc, Φn, ones(2))Finite element models are available for the 1D continuous systems defined in Section 1.1.
Data type
OneDMesh
OneDMesh(model, xmin, Nelt, bc)Construct a mesh for a beam with Nelt elements, length L and starting at xmin.
Constructor parameters
model: Structure containing the data related to the 1D system
xmin: starting position of the beam
Nelt: number of elements
bc: Boundary conditions type
:CC: Clamped - Clamped
:FF: Free - Free
:CF: Clamped - Free
:SS: Simply Supported - Simply Supported (specific to beam)
:CS: Clamped - Simply Supported (specific to beam)
:SF: Simply Supported - Free (specific to beam)
Fields
xmin: Starting position of the beam
L: Length of the beam
Nodes: Nodes of the mesh
Elt: Elements of the mesh
Ndof_per_node`: Number of degrees of freedom per node
elem_size: Size of the elements
constrained_dofs: Constrained degrees of freedom
free_dofs: Free degrees of freedom
Related functions
assemby
assembly(model::OneDstructure, mesh::OneDMesh)Compute the global stiffness and mass matrices for a 1D structure with a given mesh.
Inputs
model: OneDStructure
Beam: Beam model
WaveEquation: Bar, Rod or String model
mesh: OneDMesh
Outputs
K: global stiffness matrix
M: global mass matrix
rayleigh_damping_matrix
rayleigh_damping_matrix(K, M, α, β)
rayleigh_damping_matrix(K, M, ω1, ω2, ξ1, ξ2)Compute the Rayleigh damping matrix for a given stiffness and mass matrices
Inputs
K: Stiffness matrix
M: Mass matrix
Construction parameters
Method 1
α: Mass proportional damping coefficient
β: Stiffness proportional damping coefficient
Method 2
ω1: First natural frequency
ω2: Second natural frequency
ξ1: Damping ratio for the first natural frequency
ξ2: Damping ratio for the second natural frequency
Output
C: Rayleigh damping matrix
modal_damping_matrix
modal_damping_matrix(M, ωn, ξn, Φn)Compute the damping matrix C from modal parameters
Inputs
M: Mass matrix
ωn: Natural angular frequencies
ξn: Damping ratios
Φn: Mass-normalized mode shapes
Output
C: Damping matrix
Selection matrix
selection_matrix(mesh, selected_dofs)Compute the selection matrix for the selected dofs.
Inputs
mesh: OneDmesh
selected_dofs: Selected dofs
Output
S: Selection matrix
apply_bc - See Section 2.2.1.
eigenmode - See Section 2.2.1.
modal_matrices - See Section 2.2.1.
modal_effective_mass - See Section 2.2.1.
# Dimensions
L = 1.
d = 3e-2
# Section features
S = π*d^2/4
Iz = π*d^4/64
# Material
E = 2.1e11
ρ = 7800.
# Computation parameters
fmax = 2000.
# Initialization of the data types
beam = Beam(L, S, Iz, E, ρ)
# Mesh definition
oned_mesh = OneDMesh(beam, 0., 20, :SS)
# Construction of K and M
Kfe, Mfe = assembly(beam, oned_mesh)
# Application of the BCs
Kbc = apply_bc(Kfe, oned_mesh)
Mbc = apply_bc(Mfe, oned_mesh)
# Computation ofthe eigenmodes of the structure
ωfe, Φfe = eigenmode(Kbc, Mbc)
# Calculation of the damping matrix
Cray = rayleigh_damping_matrix(Kbc, Mbc, 1., 1.)
Cmodal = modal_damping_matrix(Mbc, ωfe, 0.01, Φfe)The state space representation of a mechanical system is expressed as: \[ \dot{\mathbf{z}}(t) = \mathbf{A}_c \mathbf{z}(t) + \mathbf{B}_c \mathbf{u}(t), \] where:
\(\mathbf{z}(t)\): State vector
\(\mathbf{u}(t)\): Input vector
\(\mathbf{A}_c\): System matrix
\(\mathbf{B}_c\): Input matrix
For a mechanical system, whose equation of motion is: \[ \mathbf{M}\ddot{\mathbf{x}}(t) + \mathbf{C}\dot{\mathbf{x}}(t) + \mathbf{K x}(t) = \mathbf{u}(t), \] the corresponding continuous-time state equation is given by: \[ \begin{bmatrix} \dot{\mathbf{x}}(t) \\ \ddot{\mathbf{x}}(t) \end{bmatrix} = \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{C} \end{bmatrix}\begin{bmatrix} \mathbf{x}(t) \\ \dot{\mathbf{x}}(t) \end{bmatrix}+ \begin{bmatrix} \mathbf{0} \\ \mathbf{M}^{-1} \end{bmatrix}\mathbf{u}(t). \]
When using a modal expansion such that \[ \mathbf{x}(t) = \mathbf{\Phi}\mathbf{q}(t), \] where \(\mathbf{\Phi}\) is the mode shapes matrix and \(\mathbf{q}(t)\) is the modal coordinate vector, a modal state space equation can be obtained. The latter is written: \[ \begin{bmatrix} \dot{\mathbf{q}}(t) \\ \ddot{\mathbf{q}}(t) \end{bmatrix} = \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{\Omega}^2 & -\mathbf{\Xi} \end{bmatrix}\begin{bmatrix} \mathbf{q}(t) \\ \dot{\mathbf{q}}(t) \end{bmatrix}+ \begin{bmatrix} \mathbf{0} \\ \mathbf{\Phi}^\mathsf{T} \end{bmatrix}\mathbf{u}(t), \] where \(\mathbf{\Omega}^2 = \text{diag}(\omega_1^2, \dots, \omega_N^2)\) and \(\mathbf{\Xi} = \text{diag}(2\xi_1\omega_1, \dots, 2\xi_N\omega_N)\).
Data type
ContinuousStateSpace
ContinuousStateSpace(Ac, Bc)Continuous-time state-space model
Fields
Ac: Continuous-time state matrix A
Bc: Continuous-time input matrix B
Related functions
ss_model
ss_model(K, M, C)Generates a continuous-time state-space model from the mass, damping, and stiffness matrices
Inputs
K: Stiffness matrix
M: Mass matrix
C: Damping matrix
Output
css: ContinuousStateSpace
ss_modal_model
ss_modal_model(ωn, ξn, ϕn)Generates a continuous-time state-space model from the mass, damping, and stiffness matrices
Inputs
ωn: Natural angular frequencies
ξn: Damping ratios
ϕn: Mass-normalized mode shapes
Output
css: ContinuousStateSpace
eigenmode
eigenmode(K, M, n = size(K, 1))Computes the eigenmodes of a system defined by its mass and stiffness matrices.
Inputs
K: Stiffness matrix
M: Mass matrix
n: Number of modes to be keep in the modal basis
Outputs
ω: Vector of natural frequencies
Φ: Mass-normalized mode shapes
modal_parameters
modal_parameters(λ)Computes the natural angular frequencies and damping ratios from the complex eigenvalues
Input
λ: Complex eigenvalues
Outputs
ωn: Natural angular frequencies
ξn: Damping ratios
c2r_modeshape
c2r_modeshape(Ψ)Converts the complex modes to real modes
Input
Ψ: Complex modes
Output
ϕn: Real modes
# System matrices
m_ss = Diagonal([2., 1.])
k_ss = [6. -2.; -2. 4.]
c_ss = [0.67 -0.11; -0.11 0.39]
# Continuous-time state space from system matrices
css = ss_model(k_ss, m_ss, c_ss)
λ, Ψ = eigenmode(css.Ac)
ω, ξ = modal_parameters(λss)
Ψr = c2r_modeshape(Ψ)
# Continuous-time state space from modal information
ωn, ϕn = eigenmode(k_ss, m_ss)
css_modal = ss_modal_model(ωn, 0.01, ϕn)Discrete-time models can be either be obtained using sampled strategies from direct-time integration methods (e.g. Newmark’s scheme). However, both approaches lead to a discrete-time state equation of the form: \[ \mathbf{z}_{k+1} = \mathbf{A}\, \mathbf{z}_k + \mathbf{B}_f\, \mathbf{u}_k + \mathbf{B}_g\, \mathbf{u}_{k+1}, \]
The previous formulation is rather non-standard, since the state vector at frequency step \(k + 1\) requires the knowledge of the input vector at time steps \(k\) and \(k + 1\). To reduce the state-space representation to its standard form, a reduced state \(\mathbf{\overline{z}}_{k+1}\) is introduced: \[ \mathbf{z}_{k+1} = \mathbf{A}\, \mathbf{x}_k + \mathbf{B}_f\, \mathbf{u}_k, \]
In doing so, the discretized state equation becomes: \[ \mathbf{\overline{z}}_{k+1} = \mathbf{A}\, \mathbf{\overline{z}}_k + \mathbf{B}\, \mathbf{u}_k \] where \(\mathbf{B} = \mathbf{B}_f + \mathbf{A\, B}_g\).
Sampling methods are based on the discretization of the solution of the continuous-time state equation. Once discretized, the state equation can be written as: \[ \mathbf{z}_{k+1} = \mathbf{A} \, \mathbf{z}_k + \mathbf{A}\, \int_{0}^{h} e^{-\mathbf{A_c}\tau}\, \mathbf{B_c}\,\mathbf{u}(\tau + kh)\, d\tau, \] where \(\mathbf{A} = e^{\mathbf{A_c} h}\) is the the discretized system matrix and \(h\) is the discretization timestep.
To cover a wide range of situations, the input vector is assumed to be parameterized as follows: \[ \mathbf{u}(\tau + kh) = f(\tau)\, \mathbf{u}_k + g(\tau)\, \mathbf{u}_{k+1}, \] where the basis function \(f(\tau)\) and \(g(\tau)\) control the evolution of the input vector between two time steps.
Introducing the parameterization of the input vector into the convolution integral, one obtains the following state equation: \[ \mathbf{z}_{k+1} = \mathbf{A}\, \mathbf{z}_k + \mathbf{B}_f\, \mathbf{u}_k + \mathbf{B}_g\, \mathbf{u}_{k+1}, \] where \(\mathbf{B}_f\) and \(\mathbf{B}_g\) are the discretized input matrices defined such that: \[ \mathbf{B}_f = \mathbf{A}\, \int_{0}^{h} e^{-\mathbf{A_c}\tau}\, \mathbf{B_c}\,f(\tau)\, d\tau \; \text{ and } \; \mathbf{B}_g = \mathbf{A}\, \int_{0}^{h} e^{-\mathbf{A_c}\tau}\, \mathbf{B_c}\,g(\tau)\, d\tau. \]
In the literature, the most commonly used sampling methods are the following1:
The Zero-Order Hold (ZOH) strategy, which assumes that the input vector is constant between two samples. This assumption is satisfied for \(f(\tau) = 1\) and \(g(\tau) = 0\) and leads to: \[ \mathbf{B} = \mathbf{B}_f = (\mathbf{A} - \mathbf{I})\, \mathbf{A}_\mathbf{c}^{-1}\, \mathbf{B_c}. \]
The First-Order Hold (FOH) sampling method, which assumes that the input vector varies linearly between two samples. This assumption is satisfied for \(f(\tau) = 1 - \tau/h\) and \(g(\tau) = \tau/h\) and leads to: \[ \begin{split} &\mathbf{B}_f + \mathbf{B}_g = (\mathbf{A} - \mathbf{I})\, \mathbf{A}_\mathbf{c}^{-1}\, \mathbf{B_c}, \\ &\mathbf{B}_g = (\mathbf{A} - \mathbf{A_c}h - \mathbf{I})\mathbf{A}_\mathbf{c}^{-2}\, \mathbf{B_c}/h. \end{split} \]
The Band-Limited Hold (BLH) sampling method, which assumes that the input signal can be approximated as a band-limited signal (i.e. the energy of the signal is concentrated in a defined frequency range). This assumption is satisfied for \(f(\tau) = h\, \delta(\tau)\) and \(g(\tau) = 0\), where \(\delta(\tau)\) is the Dirac delta function and leads to: \[ \mathbf{B} = \mathbf{B}_f = \mathbf{A}\, \mathbf{B_c}\, h. \]
Despite several methods can be found in the literature, such as the Newmark’s family schemes, this package only provides the Runge-Kutta approach, because it is an explicit method (i.e. does not require any matrix inversion). After some calculation not detailed here2, the discrete system and input matrices are expressed as: \[ \begin{split} &\mathbf{A} = \frac{1}{24}\left[24\,(\mathbf{I} + \mathbf{A_c}\, h) + 12\, \mathbf{A}_\mathbf{c}^2\, h^2 + 4\, \mathbf{A}_\mathbf{c}^3\, h^3 + \mathbf{A}_\mathbf{c}^4\, h^4\right], \\ &\mathbf{B}_f = \frac{h}{24}\left[12\, \mathbf{I} + 8\, \mathbf{A_c}\, h + 3 \mathbf{A}_\mathbf{c}^2\, h^2 + \mathbf{A}_\mathbf{c}^3\, h^3\right]\mathbf{B_c}, \\ &\mathbf{B}_g = \frac{h}{24}\left[12\, \mathbf{I} + 4\, \mathbf{A_c}\, h + \mathbf{A}_\mathbf{c}^2\, h^2\right]\mathbf{B_c}. \end{split} \]
Data type
DiscreteStateSpace
DiscreteStateSpace(Ad, Bd, Bdp)Discrete-time state-space model
Fields
Ad: Discrete-time state matrix A
Bd: Discrete-time input matrix B
Bdp: Discrete-time input matrix Bp (only for :foh method)
Related function
c2d
c2d(css::ContinuouStateSpace, h::Float64, method::Symbol)Converts a continuous-time state-space model to a discrete-time state-space model.
Inputs
css: Continuous-time state-space model
h: Sampling time.
method: Discretization method
:zoh: Zero-order Hold method
:foh: First-order Hold method
:blh: Band-limited Hold method
:rk4: 4th order Runge-Kutta method
Output
DiscreteStateSpace: Discrete-time state-space model
Ad: Discrete-time state-space matrix A
Bd: Discrete-time state-space matrix B
Bdp: Discrete-time state-space matrix Bp (only for :foh and :rk4 methods)
dss = c2d(css, 0.01, :zoh)D. Bernal. “Optimal discrete to continuous transfer for band limited inputs”, Journal of Engineering Mechanics, vol. 133 (12), pp. 1370-1377, 2007.↩︎
For further details, see: J. Ghibaudo. “Inverse estimation of sparse mechanical excitation sources by Bayesian filtering”, PhD thesis, Conservatoire national des arts et métiers, 2024.↩︎