Newmark method

\begin{equation}^{t+\Delta t} \dot U = ^t \dot U +\left[(1-\gamma)\cdot ^t\ddot U +\gamma \cdot ^{t+\Delta t} \ddot U\right]\Delta t \label{vel}\end{equation}

\begin{equation}^{t+\Delta t} U = ^t U +^t \dot U \Delta t +\left[(\frac{1}{2}-\beta)\cdot ^t\ddot U +\beta \cdot ^{t+\Delta t} \ddot U\right]\Delta t^2 \label{disp}\end{equation} rearrange Eq. \ref{disp},可以得到

\begin{equation}^{t+\Delta t} \ddot U = \frac{1}{\beta \Delta t^2}\left[ ^{t+\Delta t} U - ^t U - ^t \dot U \Delta t -^t \ddot U \Delta t^2(\frac{1}{2}-\beta)\right] \label{acc}\end{equation}

substitute Eq. \ref{acc} into Eq. \ref{vel},可以得到 \begin{equation}^{t+\Delta t} \dot U = ^{t+\Delta t}U\frac{\gamma}{\beta\Delta t} - ^{t}U\frac{\gamma}{\beta\Delta t} + ^t \dot U(1-\frac{\gamma}{\beta})+^t\ddot U (1-\frac{\gamma}{2\beta}\Delta t) \label{vell}\end{equation}

substitute Eq. \ref{acc} and Eq. \ref{vell} into the Eq. \ref{EOM} \begin{equation}M \cdot ^{t+\Delta t}\ddot U + C \cdot ^{t+\Delta t}\dot U+ K \cdot ^{t+\Delta t} U = ^{t+\Delta t} R \label{EOM}\end{equation}

可以得到Eq. \ref{disp1}或者Eq. \ref{disp2} \begin{equation}\begin{split} ^{t+\Delta t} U \left(\frac{1}{\beta\Delta t^2} M +\frac{\gamma}{\beta\Delta t} C+ K \right)=&^{t+\Delta t} R+\\ &M\left[^t U\frac{1}{\beta\Delta t^2}+^t\dot U \frac{1}{\beta\Delta t}+^t\ddot U(\frac{1}{2\beta}-1)\right]+\\ &C\left[^t U\frac{\gamma}{\beta\Delta t}+^t\dot U( \frac{\gamma}{\beta}-1)+^t\ddot U(\frac{\gamma}{2\beta}-1)\Delta t\right] \end{split} \label{disp1}\end{equation}

\begin{equation}\begin{split} ^{t+\Delta t} U \left(\frac{1}{\beta\Delta t^2} M +\frac{\gamma}{\beta\Delta t} C+ K \right)=&^{t+\Delta t} R+\\ &^t U\left(\frac{M}{\beta \Delta t^2}+\frac{C \gamma}{\beta\Delta t}\right)+\\ &^t \dot U\left(\frac{M}{\beta\Delta t}+C(\frac{\gamma}{\beta}-1)\right)+\\ &^t \ddot U\left(M(\frac{1}{2\beta}-1)+C\Delta t(\frac{\gamma}{2\beta}-1)\right) \end{split} \label{disp2}\end{equation}

\section{Exlicit method}

初值计算:已知初始位移 $^0 U$、初始速度 $^0 \dot U$和初始载荷 $^0 R$

  • 通过Eq. \ref{EOM}可以计算 $^0 \ddot U$ $$^0 \ddot U = \left(^0 F -^0 \dot U C -^0 U K\right)/M$$

由t时刻推 $t+\Delta t$时刻:已知由t时刻的位移、速度、加速度和 $t+\Delta t$时刻的载荷,

  • 通过Eq. \ref{disp1}计算 $t+\Delta t$时刻的位移;

  • 通过Eq. \ref{vell}计算 $t+\Delta t$时刻的速度;

  • 通过Eq. \ref{acc}计算 $t+\Delta t$时刻的加速度;

\section{Implicit method}

该方法基于物体在外力作用下的平衡方程,其中R为externally applied nodal point forces, F 为nodal point forces that correspond to the element stresses. F可以通过 $^t U$进行计算,在<finite element procedures in engineering analysis>section 6.3中有讨论。 \begin{equation}M ^t \ddot U = ^t R -^t F\end{equation}

时间步进方法的核心是:已知t时刻的解,求解 $t+\Delta$时刻满足Eq. \ref{gov}的解。 \begin{equation}M ^{t+\Delta t} \ddot U = ^{t+\Delta t} R - ^{t+\Delta t} F \label{gov}\end{equation} 因为t时刻的解已知,所以我们将 $^{t+\Delta t} F$表示为Eq. \ref{F},其中 $\Delta F$为由于element displacement and stress 增量引起的nodal point forces的增量。 \begin{equation}^{t+\Delta t} F =^t F + \Delta F \label{F}\end{equation}

$\Delta F$的一阶近似为Eq. \ref{app},其中 $^t K$ is a tangent stiffness matrix corresponds to the geometric and material conditions at t,见Eq. \ref{K}. $\Delta U$是nodal point displacement 的增量。 \begin{equation}\Delta F \doteq ^t K \Delta U \label{app}\end{equation}

\begin{equation}^t K =\frac{\partial ^t F}{\partial ^t U} \label{K}\end{equation} 将Eq. \ref{app}和Eq. \ref{F}带入到Eq. \ref{gov}中,得到 \begin{equation}M ^{t+\Delta t} \ddot U + ^t K \Delta U = ^{t+\Delta t} R - ^{t} F \label{gov2}\end{equation} 对于静态系统,Eq. \ref{gov2}退化为Eq. \ref{gov3}。如果系统是线性的,即Eq. \ref{K}中K独立于位移U,那么由Eq. \ref{gov3}可以直接结算 $\Delta U$,从而由 $^{t+\Delta t} U =^t U +\Delta U$得到 $t+\Delta t$时刻的解 $^{t+\Delta t}U$ \begin{equation} ^t K \Delta U = ^{t+\Delta t} R - ^{t} F \label{gov3}\end{equation} 实际问题中K与U是相关的,多余Eq. \ref{app}中的近似会造成比较大的误差,因此需要采用Newton-Raphson迭代方法获得较为准确的解。 为了获得更具有普适性的模型,这里直接对动态系统进行分析(省略了damping forces)。牛顿法的处理方式为 \begin{equation}M ^{t+\Delta t}\ddot U^{(k)} +^{t+\Delta} K^{(k-1)} \Delta U^{(k)} = ^{t+\Delta t} R -^{t+\Delta t} F^{(k-1)} \label{govd}\end{equation} 注意,对比Eq. \ref{gov2}和Eq. \ref{govd},会发现K和F左上角的t变成了 $t+\Delta t$,这是因为Eq. \ref{gov2}适用的是线性系统,在计算 $t+\Delta$时刻的U时,用的是t时刻的矩阵K,通过Eq. \ref{K}获得。而Eq. \ref{govd}针对的是非线性系统,在计算 $^{t+\Delta} U$的过程中,K不是直接由t时刻的Eq. \ref{K}获得的,因此,更适合使用上标 $^{t+\Delta}K$表示在计算 $^{t+\Delta} U$的过程中使用的K。F同理。

从第k-1次迭代到k次迭代,有Eq. \ref{U}. \begin{equation}^{t+\Delta t} U^{(k)} = ^{t+\Delta t} U^{(k-1)} + \Delta U ^{(k)} \label{U}\end{equation}

另外,根据trapezoidal time integration, 有关系式Eq. \ref{diss}和Eq. \ref{v} \begin{equation}^{t+\Delta t} U = ^t U + \frac{\Delta t}{2}\left(^t \dot U +^{t+\Delta t}\dot U \right) \label{diss}\end{equation}

\begin{equation}^{t+\Delta t} \dot U = ^t \dot U + \frac{\Delta t}{2}\left(^t \ddot U +^{t+\Delta t}\ddot U \right) \label{v}\end{equation} 结合Eq. \ref{U}, Eq. \ref{diss}和 Eq. \ref{v}可以得到Eq. \ref{ac}. \begin{equation}^{t+\Delta t}\ddot U ^{(k)} = \frac{4}{\Delta t^2}\left(^{t+\Delta t} U^{(k-1)} - ^t U +\Delta U^{(k)}\right) -\frac{4}{\Delta t} ^t U^{(k)} - ^t \ddot U \label{ac}\end{equation}

将Eq. \ref{ac}带入到Eq. \ref{govd}中,可以得到Eq. \ref{Delta}: \begin{equation}\begin{split} \left(^t K + \frac{4}{\Delta t^2} M \right)\Delta U^{(k)} =& ^{t+\Delta t} R - ^{t+\Delta t}F^{(k-1)} \\ & - M\left[\frac{4}{\Delta t^2} \left(^{t+\Delta t}U^{(k-1)} - ^t U\right)-\frac{4}{\Delta t} ^t \dot U -^t \ddot U\right] \end{split} \label{Delta}\end{equation} 算法流程如下:

  1. 根据初始条件 $^{t+\Delta t}U^{(0)} = ^t U$, $^{t+\Delta t}\dot U^{(0)} = ^t U$, $^{t+\Delta t}\ddot U^{(0)} = ^t U$, $^{t+\Delta t} K^{(0)} = ^t K$, $^{t+\Delta t}F^{(0)} = ^t F$,通过Eq. \ref{Delta}计算 $\Delta U^{(1)}$

  2. 由Eq. \ref{U} 计算 $^{t+\Delta t} U^{(1)}$

  3. $^{t+\Delta t} U^{(1)}$计算 $^{t+\Delta t} F^{(1)}$

  4. 接下来通过Eq. \ref{Delta}计算 $\Delta U^{(2)}$重复步骤231,直到 $\Delta U^{(k)}$足够小,结束循环。令 $^{t+\Delta t} U = ^{t+\Delta t} U^{(k)}$

\section{Explicit Central Difference Method}

中心差分法不涉及变量 $\beta$ $\gamma$,而是使用简单的中心差分方法: \begin{equation}^t \ddot U = \frac{1}{\Delta t^2}\left(^{t+\Delta t} U -2 ^t U + ^{t-\Delta t} U\right)\end{equation} \label{CDacc}

\begin{equation}^t \dot U = \frac{1}{2\Delta t}\left(^{t+\Delta t} U - ^{t-\Delta t} U\right)\end{equation} \label{CDvel}

将Eq. \ref{CDacc}和Eq. \ref{CDvel}带入Eq. \ref{EOM2}中,可以得到Eq. \ref{CD}。简记为Eq. \ref{CD2} \begin{equation}M \cdot ^{t}\ddot U + C \cdot ^{t}\dot U+ K \cdot ^{t} U = ^{t} R \label{EOM2}\end{equation}

\begin{equation} ^{t+\Delta t} U \left(\frac{1}{\Delta t^2} M + \frac{1}{2\Delta t}C\right) = ^{t} R - ^t U \left(K - \frac{2}{\Delta t^2}\right) - ^{t-\Delta t} U\left(\frac{1}{\Delta t^2} M -\frac{1}{2 \Delta t} C\right) \label{CD}\end{equation}

\begin{equation}^{t+\Delta t} \hat M = ^{t} \hat R\end{equation} \label{CD2}

注意初始时刻的计算需要用到 $^{0-\Delta t} U$,我们利用初始条件计算 $^{0-\Delta t} U$,见Eq. \ref{ini}。注意Central Difference Method这种方法有条件稳定, $\Delta t$不能太大,需要满足CFL稳定条件。 \begin{equation}^{0-\Delta t} U = ^0 U - ^0 \dot U \Delta t + ^0 \ddot U \frac{\Delta t^2}{2}\end{equation} \label{ini}

No comment found.

Add a comment

You must login to add a comment.

Site Maintenance

Our platform is currently undergoing maintenance. We apologize for any inconvenience. Please check back later.