Linear Mechanics

CDM(K, M, C, f, bc, u0, v0, n, Δt)
CDM(K, M, f, bc, u0, v0, n, Δt)
CDM(K, M, C, f, u0, v0, n, Δt; support=[])
CDM(K, M, f, u0, v0, n, Δt; support=[])

Solves a transient structural dynamic problem using the central difference method (CDM), an explicit time integration scheme. (support ≡ bc)

The semi-discrete system

M * ü(t) + C * u̇(t) + K * u(t) = f(t)

is advanced in time using central finite differences. The method supports:

• time-independent or time-dependent load vectors f
• time-independent or time-dependent displacement boundary conditions
• consistent treatment of constraint-induced inertia and damping terms
• vector-valued unknowns (VectorField)

Boundary conditions are applied solver-side: on constrained DOFs the prescribed displacements override the initial conditions, while on free DOFs the initial displacement u0 and velocity v0 are used. Velocities on constrained DOFs are derived from the prescribed displacements.

If f.nsteps == 1, the load is treated as time-independent. If f.nsteps == n, the load is applied time step–by–time step.

Arguments

• K::SystemMatrix Global stiffness matrix.
• M::SystemMatrix Global mass matrix (assumed diagonal or lumped for efficiency).
• C::SystemMatrix Global damping matrix. If omitted, zero damping is assumed.
• f::VectorField External load vector (time-independent or time-dependent).
• bc::Vector{BoundaryCondition} Displacement-type boundary conditions (possibly time-dependent).
• u0::VectorField Initial displacement field. Overridden on constrained DOFs.
• v0::VectorField Initial velocity field. Overridden on constrained DOFs.
• n::Int Number of time steps.
• Δt::Float64 Time step size.

Returns

• u::VectorField Displacement field at all time steps (ndof × nsteps).
• v::VectorField Velocity field at all time steps (ndof × nsteps).

The associated time vector is

t = 0 : Δt : (n-1)Δt

Notes

• The method is conditionally stable. The critical time step is governed by the highest eigenfrequency and damping:

Δt_max = T_min / π * (√(1 + ξ_max^2) - ξ_max)

where T_min is the smallest modal period and ξ_max is the largest modal damping ratio.

• The algorithm itself is agnostic to the physical meaning of the displacement field;

it operates purely on the algebraic system defined by M, C, and K.

CDMaccuracyAnalysis(ωₘᵢₙ, ωₘₐₓ, Δt, type; n=100, α=..., ξ=..., β=..., show_β=..., show_ξ=...)

Gives some functions (graphs) for accuracy analysis of the CDM method. ωₘᵢₙ and ωₘₐₓ are the square root of smallest and largest eigenvalues of the Kϕ=ω²Mϕ eigenvalue problem, Δt is the time step size. type is one of the following values:

• :SR: spectral radius
• :PDR: physical damping ratio
• :ADR: algorithmic damping ratio
• :PE: period error

For details see ^[3]. n is the number of points in the graph. The damping matrix is assembled in the following ways: C=αM+βK or C=αM+β₁K+β₂KM⁻¹K+β₃KM⁻¹KM⁻¹K+⋅⋅⋅. The latter corresponds to the damping characteristic characterized by a power series consisting of powers of the natural frequencies with odd exponents. ξᵢ (ξ in the argument list) are the values ​​of the individual members of the series corresponding to the ωₘₐₓ value. βᵢ (β in the argument list) are the coefficients of the series. (see ^[4]) Either ξ or β must be specified. ξ or β are scalars or vectors. If show_β or show_ξ is true, the corresponding β or ξ values will be sent to the output. Returns a tuple of x and y values of the graph. (Can be plotted with plot(xy))

Return: xy

Types:

• ωₘᵢₙ: Float64
• ωₘₐₓ: Float64
• Δt: Float64
• n: Int64
• α: Float64
• β: Float64 of Vector{Float64}
• ξ: Float64 of Vector{Float64}
• show_β: Boolean
• show_ξ: Boolean
• xy: Tuple{Vector{Float64},Vector{Float64}}

HHT(K, M, f, bc, u0, v0, n, Δt; α=0.0, δ=0.0, γ=0.5+δ, β=0.25*(0.5+γ)^2)
HHT(K, M, f, u0, v0, n, Δt; α=0.0, δ=0.0, γ=0.5+δ, β=0.25*(0.5+γ)^2)

Solves a transient structural dynamic problem using the Hilber–Hughes–Taylor (HHT-α) method^[1], an implicit time integration scheme with controllable numerical dissipation.

The semi-discrete system

M * ü(t) + K * u(t) = f(t)

is integrated in time using the generalized Newmark formulation with the HHT-α modification. The method supports:

• time-independent or time-dependent load vectors f
• time-independent or time-dependent displacement boundary conditions
• solver-side enforcement of Dirichlet constraints
• consistent treatment of constraint-induced inertia forces
• vector-valued unknowns (VectorField)

Boundary conditions are applied solver-side: on constrained DOFs the prescribed displacements override the initial conditions, while on free DOFs the initial displacement u0 and velocity v0 are used. Prescribed displacements are inserted at every time step, and the corresponding inertial contributions are accounted for automatically through the global system.

If f.nsteps == 1, the load is treated as time-independent. If f.nsteps == n, the load is applied time step–by–time step.

The parameters α, β, and γ control numerical dissipation and stability. If δ is given, the standard relations

γ = 0.5 + δ
β = 0.25 * (0.5 + γ)^2

are used.

Arguments

• K::SystemMatrix Global stiffness matrix.
• M::SystemMatrix Global mass matrix.
• f::VectorField External load vector (time-independent or time-dependent).
• bc::Vector{BoundaryCondition} Displacement-type boundary conditions (possibly time-dependent).
• u0::VectorField Initial displacement field. Overridden on constrained DOFs.
• v0::VectorField Initial velocity field. Overridden on constrained DOFs.
• n::Int Number of time steps.
• Δt::Float64 Time step size.

Returns

• u::VectorField Displacement field at all time steps (ndof × nsteps).
• v::VectorField Velocity field at all time steps (ndof × nsteps).

The associated time vector is

t = 0 : Δt : (n-1)Δt

Notes

• The HHT-α method is unconditionally stable for appropriate parameter choices and introduces numerical dissipation primarily in the high-frequency range.
• The algorithm itself is independent of the physical interpretation of the displacement field and operates purely on the algebraic system defined by M, K, and f.
• Damping matrices are not included explicitly; numerical dissipation is controlled via the HHT-α parameters.

Reference

Hilber, H. M., Hughes, T. J. R., Taylor, R. L. Improved numerical dissipation for time integration algorithms in structural dynamics, Earthquake Engineering & Structural Dynamics, 5(3), 283–292, 1977.

HHTaccuracyAnalysis(ωₘᵢₙ, ωₘₐₓ, Δt, type; n=100, α=0.0, δ=0.0, γ=0.5 + δ, β=0.25 * (0.5 + γ)^2)

Gives some functions (graphs) for accuracy analysis of the HHT-α method^[1]. ωₘᵢₙ and ωₘₐₓ are the square root of smallest and largest eigenvalues of the Kϕ=ω²Mϕ eigenvalue problem, Δt is the time step size. type is one of the following values:

• :SR: spectral radius
• :ADR: algorithmic damping ratio
• :PE: period error

For details see ^[2] and ^[3]. n is the number of points in the graph. For the meaning of α, β and γ see ^[1]. If δ is given, γ=0.5+δ and β=0.25⋅(0.5+γ)². Returns a tuple of x and y values of the graph. (Can be plotted with plot(xy))

Return: xy

Types:

• ωₘᵢₙ: Float64
• ωₘₐₓ: Float64
• Δt: Float64
• n: Int64
• α: Float64
• β: Float64
• γ: Float64
• δ: Float64
• xy: Tuple{Vector{Float64},Vector{Float64}}

applyBoundaryConditions!(field::Union{ScalarField,VectorField,TensorField}, supports::Vector{BoundaryCondition})

Applies boundary conditions supports on a ScalarField, VectorField or TensorField field. Mesh details are in problem. supports is BoundaryCondition.

Returns: nothing

Types:

• problem: Problem
• field: ScalarField, VectorField or TensorField
• supports: Vector{BoundaryCondition}

applyBoundaryConditions(problem::Problem, supports::Vector{BoundaryCondition})

Create a ScalarField, VectorField or TensorField and applies the boundary conditions supports on it. Mesh details are in problem. supports is BoundaryCondition.

Returns: ScalarField, VectorField or TensorField

Types:

• problem: Problem
• supports: Vector{BoundaryCondition}

applyElasticSupport!(stiffMat, elastSupp)

Applies elastic support boundary conditions elastSupp on a stiffness matrix stiffMat. Mesh details are in problem. elastSupp is a tuple of name of physical group and prescribed kx, ky and kz stiffnesses.

Returns: nothing

Types:

• stiffMat: SystemMatrix
• elastSupp: Vector{BoundaryCondition}}

dampingMatrix(K, M, ωₘₐₓ; α=0.0, ξ=..., β=...)

Generates the damping matrix for proportional damping. C = αM + βK, or C = αM + β₁K + β₂KM⁻¹K + β₃KM⁻¹KM⁻¹K + ⋯. The latter corresponds to a damping characteristic given by a power series in the natural frequencies with odd exponents. ξᵢ (ξ in the arguments) are the values of the individual terms of the series at ωₘₐₓ. βᵢ (β in the arguments) are the coefficients of the series. Either ξ or β must be specified; each may be a scalar or a vector. K is the stiffness matrix, M is the mass matrix, and ωₘₐₓ is the largest natural frequency.

Returns: dampingMatrix

Types:

• K: SystemMatrix
• M: SystemMatrix
• ωₘₐₓ: Float64
• α: Float64
• ξ: Float64 or Vector{Float64}
• β: Float64 or Vector{Float64}
• dampingMatrix: SystemMatrix

Examples

K = stiffnessMatrix(problem)
M = massMatrix(problem; lumped=true)
ωmax = 2π * 1000
C = dampingMatrix(K, M, ωmax; α=0.0, ξ=[0.02, 0.02])

elasticSupportMatrix(problem, elSupp)

Solves the elastic support matrix of the problem. elSupp is a vector of elastic supports defined in function elasticSupport. With the displacementent vector q in hand the reaction force vector fR arising from the elastic support can be solved. (fR = heatConvMat * q)

Return: elSuppMat

Types:

• problem: Problem
• elSupp: Vector{BoundaryCondition}
• elSuppMat: SystemMatrix

initialDisplacement!(u0, name; ux=..., uy=..., uz=...)

Changes the displacement values to ux, uy and uz (depending on the dimension of the problem) at nodes belonging to physical group name. Original values are in displacement vector u0.

Return: nothing

Types:

• name: String
• ux: Float64
• uy: Float64
• uz: Float64
• u0: VectorField

initialDisplacement(problem, name; ux=..., uy=..., uz=...)

Sets the displacement values ux, uy and uz (depending on the dimension of the problem) at nodes belonging to physical group name. Returns the initial displacement vector u0.

Return: u0

Types:

• problem: Problem
• name: String
• u0: VectorField
• ux: Float64
• uy: Float64
• uz: Float64

initialVelocity!(v0, name; vx=..., vy=..., vz=...)

Changes the velocity values vx, vy and vz (depending on the dimension of the problem) at nodes belonging to physical group name. Original values are in velocity vector v0.

Returns: nothing

Types:

• name: String
• v0: VectorField
• vx: Float64
• vy: Float64
• vz: Float64

initialVelocity(problem, name; vx=..., vy=..., vz=...)

Sets the velocity values vx, vy and vz (depending on the dimension of the problem) at nodes belonging to physical group name. Returns the initial velocity vector v0.

Return: v0

Types:

• problem: Problem
• name: String
• vx: Float64
• vy: Float64
• vz: Float64
• v0: VectorField

largestEigenValue(K, M)

Solves the smallest eigenvalue for a transient problem given by stiffness (heat conduction) matrix K and the mass (heat capacity) matrix M (C).

Return: λₘᵢₙ

Types:

• K: SystemMatrix
• M: SystemMatrix
• λₘᵢₙ: Float64

largestPeriodTime(K, M, bc)

Solves the largest period of time for a dynamic problem given by stiffness matrix K and the mass matrix M, bc is a vector of BoundaryCondition where the displacement is given.

Return: Δt

Types:

• K: SystemMatrix
• M: SystemMatrix
• bc: Vector{BoundaryCondition}
• Δt: Float64

loadVector(problem, loads; F=nothing)

Assembles the right-hand-side vector associated with external mechanical loads in weak-form finite element problems.

Depending on the problem dimension and the dimension of the physical group, this function handles concentrated forces, line/surface/volume loads, tractions, pressures, and body forces in 1D, 2D, and 3D.

The interpretation of the physical group dimension is:

• 1D problems:
  • Point → concentrated force
  • Line → distributed force
• 2D problems:
  • Point → concentrated force
  • Line → surface force
  • Surface → body force
• 3D problems:
  • Point → concentrated force
  • Line → edge force
  • Surface → surface force
  • Volume → body force

Load components may be specified as constants, spatial functions, or ScalarFields (nodal type). Axisymmetric and other weighted formulations can be handled by including the appropriate geometric factor (e.g. 2πr) in the load definition or coefficient.

Returns a VectorField for vector-valued problems (pdim > 1) and a ScalarField for scalar problems (pdim == 1).

If F is the deformation gradient (a TensorField), then the function solves the follower load for large displacement problems (in 3D).

Return: loadVec

Types:

• problem: Problem
• loads: Vector{BoundaryCondition}
• loadVec: VectorField
• F: Union{Nothing,TensorField}

massMatrix(problem; lumped=...)

Solves the mass matrix of the problem. If lumped is true, computes the lumped mass matrix.

Returns: massMat

Types:

• problem: Problem
• lumped: Boolean
• massMat: SystemMatrix

Examples

M = massMatrix(problem; lumped=true)

nonLinearStiffnessMatrix(problem, q)

Solves the nonlinear stiffness matrix of the problem. q is a displacement field.

Returns: stiffMat

Types:

• problem: Problem
• q: VectorField
• stiffMat: SystemMatrix

printElasticConstantsTable(; io::IO=stdout)

Prints a compact reference table of the most common relations between isotropic linear-elastic material constants:

• Young's modulus E
• Poisson's ratio ν
• Shear modulus G (also denoted as μ)
• Bulk modulus K
• Lamé's first parameter λ
• Lamé's second parameter μ (= G)

The table is meant as a quick cheat-sheet for documentation / teaching.

smallestEigenValue(K, M)

Solves the largest eigenvalue for a transient problem given by stiffness (heat conduction) matrix K and the mass (heat capacity) matrix M (C).

Return: λₘₐₓ

Types:

• K: SystemMatrix
• M: SystemMatrix
• λₘₐₓ: Float64

smallestPeriodTime(K, M, support=[])

Solves the smallest period of time for a dynamic problem given by stiffness matrix K and the mass matrix M, support is a vector of BoundaryCondition where the displacement is given.

Return: Δt

Types:

• K: SystemMatrix
• M: SystemMatrix
• support: Vector{BoundaryCondition}
• Δt: Float64

solveAxialForce(u::VectorField)

Compute axial (bar/truss) forces from a displacement field.

The input displacement field u must be nodal (VectorField), typically containing the nodal displacements of a truss or bar structure. The output is a scalar field (ScalarField), where each value represents the axial force in a truss element.

Arguments

• u::VectorField: nodal displacement field.

Returns

• ScalarField: axial forces defined per element.

Examples

u = solveDisplacement(problem, load=[loads], support=[supports])  # VectorField of nodal displacements
N = solveAxialForce(u)                               # ScalarField of axial element forces

solveBuckling(problem; load=[], support=[], n=6)

Solves the multipliers for the first n critical forces and the corresponding buckling shapes for the instability of the problem, when loads and support are applied. Result can be presented by showBucklingResults function. loads and support can be defined by load and displacementConstraint functions, respectively.

Return: buckling

Types:

• problem: Problem
• load: Vector{BoundaryCondition}
• support: Vector{BoundaryCondition}
• n: Int64
• buckling: Eigen

solveBucklingModes(K, Knl; n=6)

Solves the critical force multipliers and buckling mode shapes of a problem given by stiffness matrix K and the nonlinear stiffness matrix Knl. n is the number of buckling modes to solve. Returns the struct of critical forces and buckling modes. Results can be presented by showBucklingResults function.

Return: modes

Types:

• K: SystemMatrix
• Knl: SystemMatrix
• n: Int64
• modes: Eigen

solveDisplacement(problem;
                  load = BoundaryCondition[],
                  support = BoundaryCondition[],
                  elasticSupport = BoundaryCondition[],
                  condensed = true,
                  iterative = false,
                  reltol = sqrt(eps()),
                  maxiter = problem.non * problem.dim,
                  preconditioner = Identity(),
                  ordering = true)

Computes the displacement vector u for a given [Problem] subject to external loads, essential supports, and optional elastic supports.

The stiffness matrix is assembled internally. Depending on the value of condensed, either a reduced system (eliminating constrained degrees of freedom) or the full system with modified rows and columns is solved.

This is the high-level, user-facing displacement solver.

Arguments

• problem::Problem: Finite element problem definition (geometry, materials, DOFs, etc.).
• load::Vector{BoundaryCondition} (keyword, optional): Prescribed nodal loads.
• support::Vector{BoundaryCondition} (keyword, optional): Essential (Dirichlet) boundary conditions.
• elasticSupport::Vector{BoundaryCondition} (keyword, optional): Elastic (spring-type) supports contributing to the stiffness matrix.
• condensed::Bool (keyword): If true, solves the reduced system using only free degrees of freedom. If false, modifies the full system matrix explicitly.
• iterative::Bool (keyword): If true, uses the conjugate gradient method.
• reltol::Real (keyword): Relative tolerance for the iterative solver.
• maxiter::Int (keyword): Maximum number of iterations for the iterative solver.
• preconditioner (keyword): Preconditioner object for iterative solution (e.g. from ilu or ichol).
• ordering::Bool (keyword): If false, disables column reordering in the direct solver (lu(A, q=nothing)).

Returns

• u::VectorField: Displacement vector defined on the problem degrees of freedom.

Notes

• If problem.type == :dummy, the function returns nothing.
• When condensed = true, constrained DOFs are eliminated algebraically.
• When condensed = false, constrained rows and columns are modified explicitly.
• This function internally calls the low-level solveDisplacement(K, f; ...) routine after assembling the stiffness matrix and load vector.

solveDisplacement(K, f; support = BoundaryCondition[],
                  iterative = false,
                  reltol = sqrt(eps()),
                  maxiter = K.model.non * K.model.dim,
                  preconditioner = Identity(),
                  ordering = true)

Solves the linear system

K * u = f

for the displacement vector u, where K is a stiffness matrix and f is the load vector. Essential (Dirichlet) boundary conditions are imposed via the support argument.

This is a low-level solver operating directly on a preassembled [SystemMatrix] and load [VectorField]. It assumes that the stiffness matrix already includes all material contributions (and optional elastic supports).

If iterative = true, the reduced system is solved using the conjugate gradient method. Otherwise, a direct solver is used.

Arguments

• K::SystemMatrix: Assembled stiffness matrix associated with a [Problem].
• f::VectorField: Load vector.
• support::Vector{BoundaryCondition} (keyword, optional): Essential boundary conditions (Dirichlet-type constraints). Default is an empty vector (no supports).
• iterative::Bool (keyword): If true, uses the conjugate gradient method.
• reltol::Real (keyword): Relative tolerance for the iterative solver.
• maxiter::Int (keyword): Maximum number of iterations for the iterative solver.
• preconditioner (keyword): Preconditioner object for iterative solution (e.g. from ilu or ichol).
• ordering::Bool (keyword): If false, disables column reordering in the direct solver (lu(A, q=nothing)).

Returns

• u::VectorField: Displacement vector satisfying the prescribed boundary conditions.

Notes

• Only the unconstrained degrees of freedom are solved for; constrained values are imposed explicitly.
• This function is primarily intended for internal or advanced use. Most users should prefer the high-level solveDisplacement(problem; ...) interface.

solveEigenModes(K, M; n=6, fₘᵢₙ=1.01)

Solves the eigen frequencies and mode shapes of a problem given by stiffness matrix K and the mass matrix M. n is the number of eigenfrequencies to solve, and solves the eigenfrequencies greater than fₘᵢₙ. Returns the struct of eigenfrequencies and eigen modes. Results can be presented by showModalResults function.

Return: modes

Types:

• K: SystemMatrix
• M: SystemMatrix
• n: Int64
• fₘᵢₙ: Float64
• modes: Eigen

solveModalAnalysis(problem; support=[]; load=[], n=6, fₘᵢₙ=0.00001, directSolver=false)

Solves the first n eigenfrequencies and the corresponding mode shapes for the problem, when load and support are applied. load and contraints are optional. Result can be presented by showModalResults function. load and support can be defined by load and displacementConstraint functions, respectively. If load are given, it solves the eigenfrequencies of a prestressed structure.

Return: modes

Types:

• problem: Problem
• load: Vector{tuples}
• support: Vector{tuples}
• n: Int64
• modes: Eigen

solveStrain(q; DoFResults=false)

Solves the strain field E from displacement vector q. Strain field is given per elements, so it usually contains jumps at the boundaries of elements. Details of mesh is available in problem. If DoFResults is true, E is a matrix with nodal results. In this case showDoFResults can be used to show the results (otherwise showStrainResults or showElementResults).

Return: E

Types:

• q: VectorField
• E: TensorField

solveStress(q; T=..., T₀=..., DoFResults=false)

Solves the stress field S from displacement vector q. Stress field is given per elements, so it usually contains jumps at the boundary of elements. Details of mesh is available in problem. If DoFResults is true, S is a matrix with nodal results. In this case showDoFResults can be used to show the results (otherwise showStressResults or showElementResults). If the T temperature field (and T₀ initial temperature field if it differs from zero) is given, the function solves also the thermal stresses.

Return: S

Types:

• q: VectorField
• T: ScalarField
• T₀: ScalarField
• S: TensorField

stiffnessMatrix(problem; forceOneThread=true)

Solves the stiffness matrix of the problem.

Returns: stiffMat

Types:

• problem: Problem
• stiffMat: SystemMatrix

Examples

K = stiffnessMatrix(problem)