newmarkbeta

README

[[file:../README.org][FrameFEA - Main README]]

* Package modal performs newmark-beta time integration

* Newmark-Beta Time Integration
We have a dynamic problem of the form:

#+begin_export latex
\begin{equation} 
\label{eq:1}
\left\lbrack M \right\rbrack\left\{ \ddot{u} \right\} + \ \left\lbrack K \right\rbrack\left\{ u \right\} = \left\{ F \right\}
\end{equation}
#+end_export

We choose to use the implicit time integration method called Newmark-Beta direct time integration where:

#+begin_export latex
\begin{equation} 
\label{eq:1.1}
\left\{ \dot{u} \right\}_{n + 1} = \ \left\{ \dot{u} \right\}_{n} + \ \left\lbrack \left( 1 - \alpha \right)\left\{ \ddot{u} \right\}_{n} + \ \alpha\left\{ \ddot{u} \right\}_{n + 1} \right\rbrack\mathrm{\Delta}t
\end{equation}
#+end_export

#+begin_export latex
\begin{equation} 
\label{eq:2}
\left\{ u \right\}_{n + 1} = \ \left\{ u \right\}_{n} + \ \left\{ \dot{u} \right\}_{n}\mathrm{\Delta}t + \ \left\lbrack \left( \frac{1}{2} - \beta \right)\left\{ \ddot{u} \right\}_{n} + \ \beta\left\{ \ddot{u} \right\}_{n + 1} \right\rbrack \left( \mathrm{\Delta}t\right)^{2} 
\end{equation}
#+end_export

By rearranging \ref{eq:1.1} and \ref{eq:2} we get:

#+begin_export latex
\begin{equation} 
\label{eq:3}
\left\lbrack \widehat{K} \right\rbrack\left\{ u \right\}_{n + 1} = \left\{ \widehat{F} \right\}
\end{equation}
#+end_export

where:

#+begin_export latex
\begin{equation} 
\label{eq:3.1}
\left\lbrack \widehat{K} \right\rbrack = \ \left\lbrack K \right\rbrack + a_{0}\ \left\lbrack M \right\rbrack
\end{equation}
#+end_export

#+begin_export latex
\begin{equation} 
\label{eq:4}
\left\{ \widehat{F} \right\} = \ \left\{ F \right\}_{n + 1} + \ \left\lbrack M \right\rbrack\ \left( a_{0}\left\{ u \right\}_{n} + a_{1}\left\{ \dot{u} \right\}_{n}\  + \ a_{2}\left\{ \ddot{u} \right\}_{n} \right)
\end{equation}
#+end_export

and:

#+begin_export latex
\begin{equation} 
\label{eq:5}
a_{0} = \ \frac{1}{\beta\mathrm{\Delta}t^{2}}\ \ ;\ \ a_{1} = \ a_{0}\mathrm{\Delta}t\ \ ;\ \ a_{2} = \ \frac{1}{2\beta} - 1
\end{equation}
#+end_export

$\left{ u \right}_{0}$ and $\left{ \dot{u} \right}_{0}$ are known initial conditions, but $\left{ \ddot{u} \right}_{0}$ is not generally known, so we solve for it using:

#+begin_export latex
\begin{equation} 
\label{eq:5.1}
\left\{ \ddot{u} \right\}_{0} = \ \left\lbrack M \right\rbrack^{- 1}\left( \left\{ F \right\}_{0} - \ \left\lbrack K \right\rbrack\left\{ u \right\}_{0} \right)
\end{equation}
#+end_export

With all of the initial conditions known, we solve Eqn. \ref{eq:3} for $\left{ u \right}_{n + 1}$ and we compute $\left{ \dot{u} \right}_{n + 1}$ and $\left{ \ddot{u} \right}_{n + 1}$ from a rearrangement of Eqn. \ref{eq:2}, i.e.

#+begin_export latex
\begin{equation} 
\label{eq:5.2}
\left\{ \ddot{u} \right\}_{n + 1} = \ a_{0}\left( \left\{ u \right\}_{n + 1} - \left\{ u \right\}_{n} \right) - a_{1}\left\{ \dot{u} \right\}_{n} - \ a_{2}\left\{ \ddot{u} \right\}_{n}
\end{equation}
#+end_export

#+begin_export latex
\begin{equation} 
\label{eq:6}
\left\{ \dot{u} \right\}_{n + 1} = \ \left\{ \dot{u} \right\}_{n} + \ a_{3}\left\{ \ddot{u} \right\}_{n} + a_{4}\left\{ \ddot{u} \right\}_{n + 1}
\end{equation}
#+end_export

where:

#+begin_export latex
\begin{equation} 
\label{eq:6.1}
a_{3} = \ \left( 1 - \alpha \right)\mathrm{\Delta}t\ \ ;\ \ a_{4} = \ \alpha\mathrm{\Delta}t
\end{equation}
#+end_export

$\alpha$ and $\beta$ are chosen to control the stability of the solution. If $\alpha = \frac{1}{2}; \beta = \frac{1}{4}$ are chosen, it leads to an unconditionally stable scheme corresponding to the constant average acceleration method. If $\alpha = \frac{1}{2}; \beta = \frac{1}{6}$ then it corresponds to the linear acceleration method.

The choice of $\Delta t$ controls the accuracy. The following formula gives an estimate of $\Delta t$ that will lead to an accurate solution:

#+begin_export latex
\begin{equation} 
\label{eq:7}
\mathrm{\Delta}t = \frac{T_{\min}}{\pi}
#+end_export

Where $T_{\min}$ is the smallest period of natural vibration associated with the problem. We already have the capability of computing natural frequencies, so getting $T_{\min}$ can be done before the time integration begins.

Documentation

Index

• type Input
  • func (inData *Input) DoF() int
  • func (inData *Input) Run(calcAccel bool) (Khat, M, Q1 [][]float64, Fhat, U1, V1, A1, R1 []float64, ...)

Types

type Input

type Input struct {
	// Run Info
	Shear bool // indicates shear deformation
	Geom  bool // indicates  geometric nonlinearity

	// Node Info
	NumNodes   int           // number of Nodes
	Xyz        []femath.Vec3 // {NumNodes} X,Y,Z node coordinates (global)
	NodeRadius []float64     // {NumNodes} node size radius, for finite sizes

	// Reaction Info
	ReactedDof   []bool // {Dof} Dof's with reactions formally "r"
	UnreactedDof []bool // {Dof} Dof's without reactions formally "q" (calculated)

	// Element Info
	NumElems     int       // Number of Frame Elements
	N1, N2       []int     // {NumElems} begin and end node numbers
	Ax, Asy, Asz []float64 // {NumElems} cross section areas, incl. shear
	Jx, Iy, Iz   []float64 // {NumElems} section inertias
	ElemRoll     []float64 // {NumElems} roll of each member, radians
	E, G         []float64 // {NumElems} elastic modulus, shear modulus
	ElemDensity  []float64 // {NumElems} member densities

	L  []float64 //  node-to-node length of each element (calculated)
	Le []float64 //  effective length, accounts for node size (calculated)

	// Load Info
	StaticTol       float64   // Convergence tolerance for the static analysis when Geom == true
	AnalyzeLoadCase bool      //  Are there mechanical loads in Mechanical Loads
	MechanicalLoads []float64 // {Dof} mechanical load vectors, all load cases

	TempLoads    bool      // Are there temperature loads in the load case
	ThermalLoads []float64 // {Dof} thermal load vectors, all load cases

	EqFMech [][]float64 //  {numElems}{12} equivalent end forces from mech loads global (calculated)
	EqFTemp [][]float64 //  {numElems}{12} equivalent end forces from temp loads global (calculated)

	// Mass Info
	LumpFlag bool      // true: lumped mass matrix or false: consistent mass matrix
	NMs      []float64 // {NumNodes} node mass for each node
	NMx      []float64 // {NumNodes} node inertia about global X axis
	NMy      []float64 // {NumNodes} node inertia about global Y axis
	NMz      []float64 // {NumNodes} node inertia about global Z axis
	EMs      []float64 // {NumElems} lumped mass for each frame element

	// Time integration parameters
	Dt    float64
	Alpha float64
	Beta  float64

	// Initial Conditions
	U0 []float64 // {Dof} prescribed node displacements
	V0 []float64 // {Dof} initial nodal velocities
	A0 []float64 // {Dof} initial nodal accelerations
	// contains filtered or unexported fields
}

Input contains the necessary input data to run a modal analysis with the exception of the stiffness matrix.

func (*Input) DoF

func (inData *Input) DoF() int

DoF returns the degrees of freedom of the model described by the ModalInput structure (6*NumNodes).

func (*Input) Run

func (inData *Input) Run(calcAccel bool) (Khat, M, Q1 [][]float64, Fhat, U1, V1, A1, R1 []float64, eqErr, rmsResid float64, ok int, finalErr error)

Run performs a modal analysis of the model, and returns the dynamic stiffness and mass matricies, and the resonant mode-shapes and frequencies.