【ANSYS Workbench 学习笔记】单自由度弹簧阻尼系统理论推导与有限元验证：自由振动 | 受迫振动 | 共振 | 模态计算 | 谐响应 | 瞬态动力学

💡 简介：通过动力学基本定律建立基本动力学模型，求解自由振动与受迫振动下的系统响应，结合有限元分析软件，模拟上述工况并对比两种方式的求解结果。
💡 软件：Ansys Workbench 2020 R2
💡 模块：Modal (模态技术)、Harmonic Response (谐响应)、Transient Structural (瞬态动力学)
💡 仓库：https://gitee.com/npc-gitee/simulation 🚀
🔎 如理解有误，望不吝指正，感谢。

一、振动微分方程求解

1.1、物体自由振动的微分方程

图1所示为单自由度弹簧阻尼系统（自由振动），其中物体的质量为 $m$，坐标原点取其质心位置，向上为正，垂直位移为 $x$（若物体位于平衡位置的正方向，位移$x$为正，若物体位于平衡位置的负方向，位移$x$为负），弹簧刚度系数为 $k$，阻尼为 $c$，不计弹簧自重。

图1. 单自由度弹簧阻尼系统（自由振动）

弹簧力和阻尼力的方向始终与位移 $x$ 的方向相反：当 $x$ 为正，即方向为正方向，则 $-kx$ 和 $-c\dot{x}$ 为负，即力的方向为负方向，当 $x$ 为负，即方向为负方向，则 $-kx$ 和 $-c\dot{x}$ 为正，即力的方向为正方向；

由动力学基本定律：$F=ma$，可得：
$$m\ddot{x}=-kx-c\dot{x}$$

结合初始条件，整理得基本动力学方程：$$\{\begin{array}{l}m\ddot{x}+c\dot{x}+kx=0\\ x(0)=x_0,\dot{x}(0)={\dot{x}}_0\end{array}\text{\hspace{1em}(1-1-1)}$$

将方程两端除以 $m$，得：
$$\ddot{x}+\frac{c}{m}\dot{x}+\frac{k}{m}x=0\text{\hspace{1em}(1-1-2)}$$

式中，令 $\frac{c}{m}=2n$，$n$ 称为阻尼常量，令 $\frac{k}{m}={\omega }_0^2$，${\omega }_0$ 称为无阻尼固有频率（常量），式(1-1-2)进一步修改为：$$\ddot{x}+2n\dot{x}+{\omega }_o^{2}x=0\text{\hspace{1em}(1-1-3)}$$

为什么 ${\omega }_0$ 称为无阻尼固有频率 ？（见下文） 🚀

式(1-1-3)为二阶常系数齐次线性微分方程。

定理：设 $x(t)=x_1(t)$ 及 $x(t)=x_2(t)$ 是方程(1)的两个解，那么，对于任何常数 $C_1$，$C_2$，$x(t)=C_1{x}_1(t)+C_2{x}_2(t)$ 仍然是(1)的解。

由此定理可知，如果能找到方程(1-1-3)的两个解 $x_1(t)$ 及 $x_2(t)$，且 $x_1(t)/x_2(t)$ 不恒等于常数，那么：
$$x(t)=C_1{x}_1(t)+C_2{x}_2(t)$$
就是含有两个任意常数的解，因而就是方程(1-1-3)的通解。

当 $r$ 为常数时，指数函数 $x=e^{rt}$ 满足 $x_1(t)$ 和 $x_2(t)$，当$r$取什么数值，使 $x=e^{rt}$ 满足方程(1-1-3)。

$$\{\begin{array}{l}\dot{x}(t)=r{e}^{rt}\\ \ddot{x}(t)=r^2{e}^{rt}\end{array}$$

代入方程(1-1-3)得：
$$r^2{e}^{rt}+2nr{e}^{rt}+{\omega }_o^2{e}^{rt}=0\to (r^2+2nr+{\omega }_o^2)e^{rt}=0\to r^2+2nr+{\omega }_o^2=0$$

解得：
$$r_{1,2}=-n\pm \sqrt{n^2-{\omega }_o^2}$$

$r$ 存在以下三种情况：

• $n^2-{\omega }_o^2>0$：该情况称为大阻尼，则$$x(t)=C_1{e}^{r_{1}t}+C_2{e}^{r_{2}t}$$
• $n^2-{\omega }_o^2=0$：该情况下，阻尼等于临界阻尼($c=2\sqrt{km}$)，这时只能得到一个特解，设$x_1(t)/x_2(t)$不恒等于常数，即$x_2(t)=e^{r_{1}t}u(t)$，代入微分方程(1-1-3)，解得$u(t)=C_a+C_{b}t$，不妨选取$u(t)=t$；从而微分方程(1-1-3)的通解为：$$x(t)=C_1{e}^{r_{1}t}+C_{2}t{e}^{r_{1}t}$$
• $n^2-{\omega }_o^2<0$：该情况称为小阻尼，得到一对共轭复根 $\alpha \pm \beta i$，利用欧拉公式，得到微分方程(1-1-3)的通解为：$$x(t)=e^{\alpha t}(C_1\cos \beta t+C_2\sin \beta t)$$
  • $\alpha =-n$，$\beta =\sqrt{{\omega }_o^2-n^2}$
  • 依据余弦函数和差公式，得：
    $$x(t)=\sqrt{C_1^2+C_2^2}{e}^{-nt}\cos (\sqrt{{\omega }_o^2-n^2}t+\theta )\text{\hspace{1em}(1-1-4)}$$

基于上面的公式推导，比较好说明为什么无阻尼的固有频率为$\sqrt{k/m}$，以及有阻尼的固有频率是什么？

什么时固有频率？固有频率通常定义为系统（在没有阻尼（或阻尼很小）的情况下）自由振荡的频率。

为什么是小阻尼情形下，那临界阻尼和大阻尼情况的呢？
由下文中临界阻尼和大阻尼情形可知，结构的振动不以振荡形式存在，而是以非振荡的形式返回到平衡位置，没有振荡特征，因此，固有频率针对于小阻尼系统。既然这样是不是存在一种结构，阻尼大于结构的临界阻尼，无论负载如何，均不会发生共振？ （电流表就是利用临界阻尼特点，使指针快速到平衡位置，倒是不知道电流表会不会共振）

无阻尼是小阻尼的一种特殊情况；
三角函数的周期为 $T=\frac{2\pi }{\omega }$，则频率为 $f=\frac{\omega }{2\pi }$

因此，从小阻尼情形下分析，由方程(1-1-4)知，质量块有阻尼自由振动的频率为：$$f=\frac{1}{2\pi }\sqrt{{\omega }_o^2-n^2}=\frac{1}{2\pi }\sqrt{\frac{k}{m}-\frac{c^2}{4m^2}}=\frac{1}{2\pi }\sqrt{\frac{k}{m}}\sqrt{1-\frac{c^2}{4km}}=\frac{1}{2\pi }{\omega }_o\sqrt{1-{\zeta }^2}$$
式中，$\zeta =\frac{c}{2\sqrt{km}}$ 称为阻尼比，即阻尼除以临界阻尼；因此，有阻尼固有频率为 $\frac{1}{2\pi }{\omega }_o\sqrt{1-{\zeta }^2}$；

当无阻尼情形下，即 $c=0$，则 $\zeta =0$，因此，无阻尼固有频率为 $\frac{1}{2\pi }{\omega }_o$。

以下分情形进行讨论。

1.1.1、大阻尼情形 $n^2-{\omega }_o^2>0$

该情形下，通解为：$$x(t)=C_1{e}^{(-n+\sqrt{n^2-{\omega }_o^2})t}+C_2{e}^{(-n-\sqrt{n^2-{\omega }_o^2})t}\text{\hspace{1em}(1-1-5)}$$

当 $t=0$，$x(0)=x_0$，$\dot{x}(0)={\dot{x}}_0$，代入方程(1-1-5)及其一阶导数方程，得：
$$\{\begin{array}{l}{x}_0=C_1+C_2\\ {\dot{x}}_0=C_1(-n+\sqrt{n^2-{\omega }_o^2})+C_2(-n-\sqrt{n^2-{\omega }_o^2})\end{array}$$

解得：
$$\{\begin{array}{l}{C}_1=\frac{{\dot{x}}_0-x_0(-n-\sqrt{n^2-{\omega }_o^2})}{2\sqrt{n^2-{\omega }_o^2}}\\ \\ C_2=\frac{-\dot{x_0}+x_0(-n+\sqrt{n^2-{\omega }_o^2})}{2\sqrt{n^2-{\omega }_o^2}}\end{array}$$

因此，微分方程的解表示为：
$$x(t)=\frac{{\dot{x}}_0-x_0(-n-\sqrt{n^2-{\omega }_o^2})}{2\sqrt{n^2-{\omega }_o^2}}{e}^{(-n+\sqrt{n^2-{\omega }_o^2})t}+\frac{-\dot{x_0}+x_0(-n+\sqrt{n^2-{\omega }_o^2})}{2\sqrt{n^2-{\omega }_o^2}}{e}^{(-n-\sqrt{n^2-{\omega }_o^2})t}$$

设 $c=16N\cdot s/m$，$k=39.48N/m$，$m=1Kg$，则$n^2=(c/(2m){)}^2=64$，${\omega }_o^2=k/m=39.48$，所以$n^2-{\omega }_0^2=4.952>0$。

同时，令 $t=0$，$x(0)=1m$，$\dot{x}(0)=0m/s$，方程所对应的图形为：（曲线绘制代码见附录） 🚀

从变化趋势可知，在大阻尼情形下，质量块直接回到并停留在平衡位置。

1.1.2、临界阻尼情形 $n^2-{\omega }_o^2=0$

该情形下，通解为：
$$x(t)=C_1{e}^{-nt}+C_{2}t{e}^{-nt}\text{\hspace{1em}(1-1-6)}$$

当 $t=0$，$x(0)=x_0$，$\dot{x}(0)={\dot{x}}_0$，代入方程(1-1-6)及其一阶导数方程，得：
$$\{\begin{array}{l}{x}_0=C_1\\ {\dot{x}}_0=C_1(-n)+C_2\end{array}$$

解得：
$$\{\begin{array}{l}{C}_1=x_0\\ C_2={\dot{x}}_0+n{x}_0\end{array}$$

因此，微分方程的解表示为：
$$x(t)=x_0{e}^{-nt}+(\dot{x}+n{x}_0)t{e}^{-nt}$$

设 $c=12.567N\cdot s/m$，$k=39.48N/m$，$m=1Kg$，则$n^2=(c/(2m){)}^2=39.48$，${\omega }_o^2=k/m=39.48$，所以$n^2-{\omega }_o^2=0$

同时，令 $t=0$，$x(0)=1m$，$\dot{x}(0)=0m/s$，方程所对应的图形为：（曲线绘制代码见附录） 🚀

在相同的初始条件下，图中为临界阻尼和大阻尼情形下的系统响应，其中虚线为大阻尼情形下系统响应，实线为临界阻尼情形下系统响应。两种变化趋势较为相似，均表现为质量块直接回到并停留在平衡位置，但是，临界阻尼回到平衡位置时间最短。

1.1.3、小阻尼情形 $n^2-{\omega }_o^2<0$

该情形下，通解为：
$$x(t)=e^{-nt}(C_1\cos (\sqrt{{\omega }_o^2-n^2}t)+C_2\sin (\sqrt{{\omega }_o^2-n^2}t))\text{\hspace{1em}(1-1-7)}$$

当 $t=0$，$x(0)=x_0$，$\dot{x}(0)={\dot{x}}_0$，代入方程(1-1-7)及其一阶导数方程，得：
$$\{\begin{array}{l}{x}_0=C_1\\ {\dot{x}}_0=-n{C}_1+C_2\sqrt{{\omega }_o^2-n^2}\end{array}$$

解得：
$$\{\begin{array}{l}{C}_1=x_0\\ \\ C_2=\frac{{\dot{x}}_0+n{x}_0}{\sqrt{{\omega }_o^2-n^2}}\end{array}$$

因此，微分方程的解表示为：
$$x(t)=e^{-nt}(x_0\cos (\sqrt{{\omega }_o^2-n^2}t)+\frac{{\dot{x}}_0+n{x}_0}{\sqrt{{\omega }_o^2-n^2}}\sin (\sqrt{{\omega }_o^2-n^2}t))$$

设 $c=2N\cdot s/m$，$k=39.48N/m$，$m=1Kg$，则$n^2=(c/(2m){)}^2=1$，${\omega }_o^2=k/m=39.48$，所以$n^2-{\omega }_o^2=-38.48<0$。

同时，令 $t=0$，$x(0)=1m$，$\dot{x}(0)=0m/s$，方程所对应的图形为：（曲线绘制代码见附录） 🚀

从变化趋势可知，在小阻尼情形下，质量块围绕平衡位置振动，振幅衰减。

1.1.4、无阻尼情形 $n=0$

$n=0$，即 $c=0$，则基本动力学方程可写成：
$$\{\begin{array}{l}m\ddot{x}+kx=0\\ x(0)=x_0,\dot{x}(0)={\dot{x}}_0\end{array}\text{\hspace{1em}(1-1-8)}$$

微分方程(1-1-8)求解可在小阻尼的基础上将 $n=0$，得：
$$x(t)=C_1\cos ({\omega }_{o}t)+C_2\sin ({\omega }_{o}t)\text{\hspace{1em}(1-1-9)}$$

当 $t=0$，$x(0)=x_0$，$\dot{x}(0)={\dot{x}}_0$，代入方程(1-1-9)及其一阶导数方程，得：
$$\{\begin{array}{l}{x}_0=C_1\\ {\dot{x}}_0=C_2{\omega }_o\end{array}$$

解得：
$$\{\begin{array}{l}{C}_1=x_0\\ \\ C_2=\frac{\dot{x_0}}{{\omega }_o}\end{array}$$

因此，微分方程的解表示为：
$$x(t)=x_0\cos ({\omega }_{o}t)+\frac{{\dot{x}}_0}{{\omega }_o}\sin ({\omega }_{o}t)$$

设 $k=39.48N/m$，$m=1Kg$，同时，令 $x(0)=1m$，$\dot{x}(0)=0m/s$，方程所对应的图形为：（曲线绘制代码见附录） 🚀

阻尼是一种能量耗散的机制，由于振动系统中阻尼为0，且没有外部能量输入，所以振动系统的运动位移没有衰减也没有增加。

1.2、物体受迫振动的微分方程

图2所示为单自由度弹簧阻尼系统（受迫振动），其中物体的质量为 $m$，坐标原点取其质心位置，向上为正，垂直位移为 $x$，弹簧刚度系数为 $k$，阻尼为 $c$，不计弹簧自重。

图2. 单自由度弹簧阻尼系统（受迫振动）

弹簧力和阻尼力始终与位移 $x$ 方向相反，激励 $F$ 的方向随时间改变。

由动力学基本定律：$F=ma$，可得：
$$m\ddot{x}=F-kx-c\dot{x}$$

结合初始条件，整理得基本动力学方程：$$\{\begin{array}{l}m\ddot{x}+c\dot{x}+kx=F_0\sin \omega t\\ x(0)=0,\dot{x}(0)=\dot{x}\end{array}\text{\hspace{1em}(1-2-1)}$$

对比自由振动和受迫振动，二者差别不大，多了外部激励，即 $F_0\sin \omega t$。

方程(1-2-1)中，$m\ddot{x}+c\dot{x}+kx=F_0\sin \omega t$ （①）称为二阶常系数非齐次线性微分方程的一般形式，$m\ddot{x}+c\dot{x}+kx=0$ （②）称为非齐次方程（①）所对应的齐次方程。

1.2.1、求解微分方程通解 $\bar{x}(t)$

定理：设$x=x^{\ast }(t)$ 是方程（①）的解，$x=\bar{x}(t)$ 是方程（②）的解，那么$$x=\bar{x}(t)+x^{\ast }(t)$$仍是方程（①）的解。

由此定理可知，需要求解$\bar{x}(t)$和$x^{\ast }(t)$，$\bar{x}(t)$是方程(1-2-1)齐次方程的通解，并且这里考虑小阻尼情形，所以可得：
$$\bar{x}(t)=e^{-nt}(C_1\cos (\sqrt{{\omega }_o^2-n^2}t)+C_2\sin (\sqrt{{\omega }_o^2-n^2}t))\text{\hspace{1em}(1-2-2)}$$

阻尼比 $\zeta =\frac{c}{2\sqrt{km}}$，阻尼常量 $n=\frac{c}{2m}$，所以$\zeta {\omega }_o=\frac{c}{2\sqrt{km}}\cdot \sqrt{\frac{k}{m}}=\frac{c}{2m}=n$，$\sqrt{{\omega }_o^2-n^2}$ 为有阻尼固有频率，记作 ${\omega }_d$，方程(1-2-2)整理得：
$$\bar{x}(t)=e^{-\zeta {\omega }_{o}t}(C_1\cos ({\omega }_{d}t)+C_2\sin ({\omega }_{d}t))\text{\hspace{1em}(1-2-3)}$$

1.2.2、求解微分方程特解 $x^{\ast }(t)$

对于形如 $f(x)=e^{\lambda x}(A\cos {\omega }_{r}t+B\sin {\omega }_{r}t)$，方程具有形如 $x^{\ast }(t)=x^k{e}^{\lambda x}(a\cos {\omega }_{r}x+b\sin {\omega }_{r}t)$的特解，而本方程中$f(x)=F_0\sin (\omega t)$符合函数形式要求，其中$\lambda =0$，$A=0$，$B=F_0$，${\omega }_r=\omega$；

$k$ 按 $\lambda +i{\omega }_r$ 按不是特征方程的根、或是特征方程的根取0或1。

基本动力学方程的特征方程为 $m{r}^2+cr+k=0$，由于 $\lambda +i{\omega }_r$ 不是特征根（等式左侧有虚部，等式右侧无虚部），所以应取 $k=0$，故设特解为：$$x^{\ast }(t)=a\cos \omega t+b\sin \omega t$$

求导得：
$$\{\begin{array}{l}{x}^{\ast }(t{)}^{{}^{\prime }}=-a\omega \sin \omega t+b\omega \cos \omega t\\ x^{\ast }(t{)}^{{}^{\prime \prime }}=-a{\omega }^2\cos \omega t-b{\omega }^2\sin \omega t\end{array}$$

代入方程(1-2-1)，整理得：
$$(-m{\omega }^{2}a+c\omega b+ka)\cos \omega t+(-m{\omega }^{2}b-c\omega a+kb)\sin \omega t=F_0\sin \omega t$$

比较两端同类项系数得：
$$\{\begin{array}{l}-m{\omega }^{2}a+c\omega b+ka=0\\ -m{\omega }^{2}b-c\omega a+kb=F_0\end{array}$$

由此解的：
$$\{\begin{array}{l}a=-\frac{F_{0}c\omega }{(k-m{\omega }^2{)}^2+c^2{\omega }^2}\\ \\ b=\frac{F_0(k-m{\omega }^2)}{(k-m{\omega }^2{)}^2+c^2{\omega }^2}\end{array}\text{\hspace{1em}(1-2-4)}$$

由前面知道 $\zeta {\omega }_0=\frac{c}{2m}$，则 $c=2\zeta {\omega }_{o}m$，所以 $c^2{\omega }^2/k^2=4{\zeta }^2{\omega }_o^2{m}^2{\omega }^2/k^2=4{\zeta }^2{\omega }_o^2{\omega }^2/{\omega }_o^4=(2\zeta s{)}^2$，其中$s=\omega /{\omega }_o$；因此，将方程(1-2-4)分母提取 $k^2$，整理得：
$$\{\begin{array}{l}a=-\frac{(F_0/k)2\zeta s}{(1-s^2{)}^2+(2\zeta s{)}^2}\\ \\ b=\frac{(F_0/k)(1-s^2)}{(1-s^2{)}^2+(2\zeta s{)}^2}\end{array}$$

因此，方程为：
$$x^{\ast }(t)=\frac{(F_0/k)(1-s^2)}{(1-s^2{)}^2+(2\zeta s{)}^2}\sin \omega t-\frac{(F_0/k)2\zeta s}{(1-s^2{)}^2+(2\zeta s{)}^2}\cos \omega t\text{\hspace{1em}(1-2-5)}$$

对于方程 $x_p(t)=C_1\sin \omega t+C_2\cos \omega t$，利用正弦函数和差公式可得：
$$x_p(t)=\sqrt{C_1^2+C_2^2}\sin (\omega t+{\theta }_o)$$

• 当 $C_1>0$，则 ${\theta }_o=\arctan \frac{C_2}{C_1}$；
• 当 $C_1<0$，$C_2>0$， 则 ${\theta }_o=\arctan \frac{C_2}{C_1}+\pi$；
• 当 $C_1<0$，$C_2<0$，则 ${\theta }_o=\arctan \frac{C_2}{C_1}-\pi$；

经计算：

• $\sqrt{a^2+b^2}=\frac{F_0/k}{\sqrt{(1-s^2{)}^2+(2\zeta s{)}^2}}=\frac{F_0}{k}\beta$，其中$\beta =\frac{1}{\sqrt{(1-s^2{)}^2+(2\zeta s{)}^2}}$ 称为振幅放大因子；
• 当$b>0$，则 ${\theta }_o=\arctan \frac{a}{b}=-\arctan \frac{2\zeta s}{1-s^2}$，记 ${\theta }_o=-\theta$；
• 当 $b<0$，则 ${\theta }_o=\arctan \frac{a}{b}-\pi =-\arctan \frac{2\zeta s}{1-s^2}-\pi$，记 ${\theta }_o=-\theta$；（$a$肯定小于0）

方程(1-2-5)整理得：
$$x^{\ast }(t)=\frac{F_0}{k}\beta \sin (\omega t-\theta )\text{\hspace{1em}(1-2-6)}$$

1.2.3、利用初始条件，求解系数 $C_1$ 和 $C_2$

联立方程(1-2-3)和方程(1-2-6)得方程(1-2-1)的通解：
$$x(t)=\bar{x}(t)+x^{\ast }(t)=e^{-\zeta {\omega }_{o}t}(C_1\cos ({\omega }_{d}t)+C_2\sin ({\omega }_{d}t))+\frac{F_0}{k}\beta \sin (\omega t-\theta )\text{\hspace{1em}(1-2-7)}$$

将方程(1-2-1)可知，$x(0)=x_0$，$\dot{x}(0)={\dot{x}}_0$，代入方程(1-2-7)得：

$$x_0=C_1-\frac{F_0}{k}\beta \sin \theta \to C_1=x_0+\frac{F_0}{k}\beta \sin \theta$$

$$\dot{x}(t)=(-\zeta {\omega }_o)e^{-\zeta {\omega }_{o}t}(C_1\cos {\omega }_{d}t+C_2\sin {\omega }_{d}t)+e^{-\zeta {\omega }_{o}t}(-C_1{\omega }_d\sin {\omega }_{d}t+C_2{\omega }_d\cos {\omega }_{d}t)+\frac{F_0}{k}\beta \omega \cos (\omega t-\theta )$$

$${\dot{x}}_0=-\zeta {\omega }_o{C}_1+C_2{\omega }_d+\frac{F_0}{k}\beta \omega \cos \theta$$

代入 $C_1$，得：
$${\dot{x}}_0=-\zeta {\omega }_o(x_0+\frac{F_0}{k}\beta \sin \theta )+C_2{\omega }_d+\frac{F_0}{k}\beta \omega \cos \theta$$

解得 $C_2$：
$$C_2=\frac{{\dot{x}}_0}{{\omega }_d}+\zeta \frac{{\omega }_o}{{\omega }_d}{x}_0+\frac{F_0}{k}\beta (\zeta \frac{{\omega }_o}{{\omega }_d}\sin \theta -\frac{\omega }{{\omega }_d}\cos \theta )$$

将 $C_1$ 和 $C_2$ 代入方程(1-2-7)，得：

$$x(t)=\stackrel{初始响应条件}{\overbrace{e^{-\zeta {\omega }_{o}t}(x_0\cos {\omega }_{d}t+\frac{{\dot{x}}_0+\zeta {\omega }_o{x}_0}{{\omega }_d}\sin {\omega }_{d}t)}}+\stackrel{自由伴随振动}{\overbrace{\frac{F_0}{k}\beta e^{-\zeta {\omega }_{o}t}(\sin \theta \cos {\omega }_{d}t+\frac{{\omega }_o}{{\omega }_d}(\zeta \sin \theta -s\cos \theta )\sin {\omega }_{d}t)}}+\stackrel{受迫振动响应}{\overbrace{\frac{F_0}{k}\beta \sin (\omega t-\theta )}}$$

其中，（虽然前面有提到，这里汇总一下）

• ${\omega }_o=\sqrt{k/m}$ 为无阻尼的固有频率；
• $s=\omega /{\omega }_o$；
• $\zeta =c/(2\sqrt{km})$ 为阻尼比；
• ${\omega }_d={\omega }_o\sqrt{1-{\zeta }^2}$ 为带阻尼的固有频率；
• $\beta =\frac{1}{\sqrt{(1-s^2{)}^2+(2\zeta s{)}^2}}$ 为振幅放大因子；
• $\theta$ 为相位角；

绘制全响应曲线前需要明确，以某种参数单位输入，输出的$x(t)$是什么单位？m? mm?

方程(1-2-1)：$m\ddot{x}+c\dot{x}+kx=F_0\sin \omega t$

• $m$ 的单位为 $Kg$；
• $c$ 的单位为 $N\cdot s/m$
• $k$ 的单位为 $N/m$；
• $F_0$ 的单位为 $N$；
• 假设 $x$ 的单位为 $m$；
• 假设 $\dot{x}$ 的单位为 $m/s$；
• 假设 $\ddot{x}$ 的单位为 $m/s^2$；

将这些单位代入方程中：$Kg\cdot \frac{m}{s^2}+\frac{N\cdot s}{m}\cdot \frac{m}{s}+\frac{N}{m}\cdot m=N+N+N=N$；符合左右两边单位相同要求。
所以，$k$、$c$、$m$（质量）以上述单位输入，输出单位为 $m$（米）。

全响应 $x(t)$ 🚀 由三部分组成，第一部分为初始响应条件，第二部分为自由伴随振动，第三部分为受迫振动响应，前两项主要由两部分乘积组成，分别为指数函数 $e^{-\zeta {\omega }_{o}t}$ （单调递减）和三角函数（在-1~1之间变化），因此随着时间的推移，这两项振动将逐渐消失，最后形成稳定的正弦曲线形式的受迫振动（方程最后一项），如图3所示。

谐响应正是研究时间与位移关系中后面段的稳定受迫振动。

绘制全响应曲线：（曲线绘制代码见附录） 🚀

图3. 各项响应曲线

图4. 全响应曲线

图3中，$x(t)-1$ 为初始响应条件；$x(t)-2$ 为自由伴随振动；$x(t)-3$ 为为受迫振动响应，其中初始响应条件与自由振动下的小阻尼振动响应 🚀是相同。

图4中，$x(t)$ 为全响应变化曲线。

1.3、 为什么激励频率与结构固有频率相近会发生振幅显著增加（共振）？

共振是一个物理学概念，指的是在特定条件下，一个系统在外部力的作用下，其振动幅度可能显著增加的现象。

随着时间的推移前两项逐渐消失，所以共振分析主要集中在受迫振动响应上：

$$x^{\ast }(t)=\frac{F_0}{k}\beta \sin (\omega t-\theta )=\frac{F_0}{k}\cdot \frac{1}{\sqrt{(1-\frac{{\omega }^2}{{\omega }_o^2}{)}^2+(2\zeta \frac{\omega }{{\omega }_o}{)}^2}}\cdot \sin (\omega t-\theta )$$

当系统确定后，其 $k$、$c$、$m$ 将保持不变，进而固有频率${\omega }_o$ 和阻尼比 $\zeta$ 也保持不变，与外界载荷无关，因此固有频率和阻尼比是系统的固有属性。对于方程而言，当外载荷以什么形式输入，$x^{\ast }(t)$振幅显著增加？

由于 $k$、${\omega }_o$、$\zeta$ 均为常量，

• 当施加较大的$F_0$，则 $x^{\ast }(t)$ 将变大；
• 当系统的阻尼比较小时候，外部载荷 $\omega$ 和固有频率 ${\omega }_o$ 相等，此时 $\sqrt{(1-\frac{{\omega }^2}{{\omega }_o^2}{)}^2+(2\zeta \frac{\omega }{{\omega }_o^2})}$ 趋近于0，那么 $x^{\ast }(t)$ 将趋于无穷大，这也解释了为什么外部激励频率与固有频率相近导致振幅显著增加。

频率不是发生共振的唯一条件，激励方向也是一个重要的因素。

二、有限元模态计算验证

模态计算是用来计算线性结构的动力学特性。

动力学特性：

• 结构的固有频率；
• 结构的模态振型；(模态分析得到的变形量不是真实变形，而是比例值，谐响应和瞬态动力学是真实变形量)
• 振型参与系数；
• 有效质量；

这里通过模态计算获得结构的固有频率与前文计算所得进行比较。

由前文叙述，无阻尼的固有频率为 ${\omega }_o=\sqrt{\frac{k}{m}}$，有阻尼的固有频率为 ${\omega }_d={\omega }_o\sqrt{1-{\zeta }^2}$。

令密度为 $1000Kg/m^3$，体积为 $0.001m^3$，则质量为 $1Kg$；弹簧系数为 $39.48N/m$，阻尼为 $2N\cdot s/m$。

计算所得：

• 无阻尼固有频率： $f_o=\frac{1}{2\pi }{\omega }_o=\frac{1}{2\pi }\sqrt{\frac{39.48N/m}{1Kg}}=1Hz$
• 有阻尼固有频率：$f_d=\frac{1}{2\pi }{\omega }_o\sqrt{1-{\zeta }^2}=\frac{1}{2\pi }\sqrt{\frac{39.48N/m}{1Kg}}\sqrt{1-(\frac{2N\cdot s/m}{2\sqrt{39.48N/m\times 1Kg}}{)}^2}=0.987Hz$

选用模块：Modal
参考视频：例10-2 二自由度弹簧振子模态分析 🚀

2.1、无阻尼

模态计算主要步骤如下：

① 双击模块【Modal】；

② 双击Engineering Data，将密度修改为$1000Kg/m^3$；

③ 右键Geometry => Replace Geometry => 选择对应模型，将三维模型添加到模块下；

④ 双击Model，进入Mechanical界面，点击Solid选择材料，选择完成后物体的质量为1Kg，设置为Rigid （刚体）；

⑤ Model => 右键Inset => Connections；

⑥ 选择Connection => 右键Inset => Spring =>

• 设置 Longitudinal Striffness 为 39.48 N/m（注意单位）；
• 设置 Scope 为 Body-Ground；
• 设置弹簧接地端的坐标，目的使弹簧朝向与连接的物体垂直；
• 设置 Face 为物体的某一个面，之后仅这个面的垂直方向平移自由度为 free，其它5个自由度均为约束；
• 设置 Behavior 为 Rigid

⑦ 选择Connection => 右键Inset => Joint => 设置约束；（另一种约束方式：见下文 🚀）

⑧ 点击 Mesh => 右键 Generate Mesh；

⑨ 点击 Modal => 右键 Solve；

⑩ 计算完成之后，点击 Solution 查看固有频率；

有限元计算得到的模态频率与公式计算结果一致。

由于只有一个自由度，所以只有一阶固有频率。

当 Joint => Type中没有 General 选项，可以选用该种约束：

点击 Modal（A5）=> 右键 Insert => Remote Displacement =>

• 选择质量块的一个面；
• 除了Y方向平移为Free，其余均为0；
• Behavior：Rigid；
  

2.2、有阻尼

① 在无阻尼的设置下，输入阻尼 $2N\cdot s/m$；

② 在Analysis Settings => Solver Controls => Damped 选择 Yes；

③ 点击Modal => 右键 Solve，待计算完成，点击Solution查看固有频率 ；

有限元计算得到的模态频率与公式计算结果一致。

PS：没设置阻尼的时候，弹簧还是弹簧形状，设置阻尼后，弹簧就变成了圆柱。

三、有限元谐响应验证

什么是谐响应分析？ 确定一个结构在已知频率的正弦（简谐）载荷作用下结构响应的技术。

如果载荷为非正弦类的周期载荷，谐响应做不了，需要在瞬态动力学模块中去完成。

激励输入：

• 已知大小和频率的谐波载荷（力、压力和强迫位移）；
• 可以输入多种类型的载荷，作用位置、先后顺序（相位）、作用类型等均可不同，但是激励频率必需相等。

谐响应载荷输入参数要求：

• 载荷的幅值 $F_{max}$；
• 载荷的相位角 $\theta$；
• 相位角表示两个或更多载荷之间的相位移动；
  • 如果只定义一个载荷，不需要定义相位角，即采用默认值0。
• 载荷的激励频率范围 $f_{min}<f<f_{max}$；（参考模态分析前几阶的频率，确定激励频率范围）

则输入激励为：
$$F_i=F_{max}\sin (2\pi ft+\theta )$$

响应输出：

• 结构输出的响应的频率和激励频率相同，比如输入载荷激励频率为1Hz，输出响应的频率也是1Hz。
  • 正如前文所述，该响应为受迫振动响应，不含初始响应条件和自由伴随振动。 由 全响应 $x(t)$ 🚀 可知，受迫振动响应的频率正好就是激励频率。
• 对于任何结构阻尼，共振点相位角为±90°。

谐响应的输出可以得到幅频曲线图（关于振幅 $z$ 和激励频率点的图）和相频曲线图（关于相位 $\psi$ 和激励频率点的图）。则某个激励下输出响应为：$$x(t)=z\sin (2\pi ft+\psi )$$

ANSYS Workbench 提供了完全法和模态叠加法（振型叠加法）两种求解方法：

• 速度较慢，求解激励频率点间隔均匀，精度不依赖于模态分析中模态求解阶数；
• 速度较快，支持共振点附近激励频率点非均匀分布求解，精度依赖于模态分析中模态求解阶数（可通过有效质量判断模态阶数是否充足）；

本实例中两种方法均进行了尝试，大差不差，毕竟只有一个自由度，也就只有一阶固有频率。

3.1、基于完全法的谐响应计算

主要步骤如下：

① 模块：Harmonic Response，这个模块既可以用完全法去计算，也可以用模态叠加法去求解。

② 材料属性：同模态计算一样，将密度设置为 $1000Kg/m^3$。

③ 模型导入：同模态计算一样，导入模型，进入Mechanical界面，选择材料，谐响应需要施加激励，而刚体无法施加，所以这里设置为柔性体。

④ 连接：同模态计算设置一样，添加 Spring (设置弹簧系数 $39.48N/m$ 和阻尼 $2N\cdot s/m$) 和 Joint。

⑤ 网格划分：为了尽可能模拟刚体效果，这里在网格划分的时候，将所有的边设置的段数为1，这样最后只有一个单元，然后一键Mesh。

⑥ 分析设置：设置激励范围、间隔点、求解方式、结构阻尼常量系数。

• Range Minimum：激励范围最小值；
• Range Maximum：激励范围最大值；
• Solution Intervals：间隔点，比如激励范围10Hz，间隔10，就是在这个频率范围内，均匀取10个激励频率点；
• Solution Method：Full，完全法；
• Constant Structural Damping Coefficient：结构阻尼常量系数$C$；
  • $C=2\zeta {\omega }_{o}m$，其中 $\zeta$ 为阻尼比，${\omega }_o$ 为固有频率，$m$ 为质量，计算所得 $C=2\times 0.16\times 1\times 1=0.32$
  • 前面弹簧处设置了阻尼，这里也设置了阻尼，系统怎么处理？
    • 如果只是设置1处，均能得到正确结果；
    • 如果两处均设置，且输入的阻尼正确，则结果错误，振幅更小，有点两处阻尼叠加效果；

⑦ 添加载荷：输入一个Force，设置Force的作用面以及力的大小，之后就可以开始计算。

⑧ 后处理：查看Bode图和振动云图。

• Solution => Insert => Frequency Response => Stress/Strain/Deformation… 查看应力/应变/变形…的幅频曲线和相频曲线（幅频图+相频图=Bode图）
• Solution => Insert => Deformation => Total，设置查看的频率，查看模型的变形云图。

图中所示为不同激励频率下的响应，从中可以知道，在1Hz的时候振幅最大，即在1Hz的激励频率下，结构的振幅存在显著增加，结构发生共振，共振的响应为：（$2\pi \ast 1t=6.283t$，$89.993°$约等于$\frac{\pi }{2}=1.571$）
$$x(t)=79.154\sin (6.283t-1.571)$$

将参数代入 全响应 $x(t)$ 🚀 中的受迫振动响应中，得到：
$$x(t)=79.578\sin (6.283t-1.571)$$

有限元计算得到的受迫振动响应与公式计算结果一致。

图5. 谐响应分析（受迫振动响应）

3.2、基于模态叠加法的谐响应计算

模态叠加法同完全法的步骤基本一致，主要步骤如下：

① 模块：Modal 和 Harmonic Response，建立二者数据关联。

② 材料属性：同完全法一样。

③ 模型导入：同完全法一样。

④ 连接：同完全法一样。

⑤ 网格划分：同完全法一样。

⑥ 模态计算：直接求解模态。

⑦ 分析设置：基本操作同完全法一样，除了：

• Solution Method，系统自动设置为 Mode Superposition 。
• Cluster Number：采用共振点附近激励频率点非均匀分布求解；

⑧ 添加载荷：同完全法一样。

⑨ 后处理：同完全法一样。

Modal 模块下的 Analysis Setting => Solver Controls => Damped：Yes 会报错：MSUP Harmonic and MSUP Transient analysis doesn't support Full Damped or Unsymmetric solver type selection of the upstream Modal analysis.

所以选No.

同完全法不同的是，在弹簧处设置阻尼，分析设置中也设置阻尼，无论设置一处还是设置两处，计算结果一样。

图中所示为不同激励频率下的响应，从中可以知道，在1Hz的时候振幅较大（最大的在其附近，可能是模态叠加法误差导致），即在1Hz的激励频率下，结构的振幅存在显著增加，结构发生共振，共振的响应为：（$2\pi \ast 1t=6.283t$，$90°$等于$\frac{\pi }{2}=1.571$）
$$x(t)=79.154\sin (6.283t-1.571)$$

同完全法的计算结果一致，同理论解同样一致。

四、有限元瞬态动力学验证

什么是瞬态动力学分析？ 它是确定随时间变化载荷作用下结构响应的技术，全过程模拟结构在初始条件和外界载荷作用下的响应，不同于谐响应，瞬态动力学的响应是全响应。

瞬态动力学同样提供了两种方法：模态叠加法和完全法，基本的操作类似，参考视频中采用模态叠加法，由于模态叠加法计算的初始条件只能是速度和位移都为0，所以本示例采用完全法。

参考视频：例11-2 二自由度弹簧振子谐响应分析 🚀

主要步骤如下：

① 模块：Transient Structural。

② 材料属性：同谐响应计算一样。

③ 模型导入：同谐响应计算一样。

④ 连接：同谐响应计算一样。

⑤ 网格划分：同谐响应计算一样。

⑥ 创建 Named Selections：选择一个Node => 右键 Create Name Selection… => 命名 => OK。

⑦ 载荷步设置：为了设置初始条件非零，需要分两个载荷步，

• 第一个载荷步：0s~0.001s，添加节点位移（Nodal Displacement），本示例设置为1m，在【Tabular Data 】中，点击 0.001s~11s 栏 => 右键 Activate/Deactivate at this step，则该栏变成灰色，表示该位移设置在第二个载荷步不起作用；
• 第二个载荷步：0.001s~11s，添加外力载荷 $\sin (360\ast (time-0.001))$（Y方向），在【Tabular Data 】中，点击 0s~0.001s 栏 => 右键 Activate/Deactivate at this step，则该栏变成灰色，表示该载荷在第一个载荷步不起作用；

初始条件设置：

从图中黄色框框中，将0s中Y方向的位移设置为0，那么从0s到0.001s的位移为1000mm，所需时间为0.001s，所以在0.001s时刻的速度为1000m/s。最后初始条件就是：位移1m，速度为1000m/s。

若将0s中Y方向的位移设置为1000mm，那么从0s到0.001s的位移为0mm，所以在0.001s时刻的速度为0m/s。最后初始条件就是：位移1m，速度为0m/s。

外力载荷设置：

外力载荷大小为1N，频率为1Hz，考虑只有一个外力，所以相位为0，加上第一个载荷步0.001s，因此为 $\sin (2\pi (time-0.001))$，由于软件中对于三角函数的计算采用的不是弧度，所以需要变成角度，即 $\sin (360\ast (time-0.001))$。

弹簧一端连着质量块，一端连着地，其中连接地的点是固定不动的定点，当弹簧的长度小于振幅，质量块将会沿着定点左右运动，这个情形会导致计算存在误差，所以需要将弹簧长度设置的最大振幅大。

弹簧长度设置：修改连接地处的坐标。

⑧ 分析设置：设置载荷步结束时间、初始时间步、最小时间步、最大时间步、时间积分。注意，第一个载荷步 Time Integration 设置成 off，第二个载荷步 Time Integration 设置成on。

⑨ 计算求解与后处理：同谐响应一样。

软件中对于系统的响应在幅值上不区分正负，因此这里将数据导出做一下处理，由于理论解和仿真解高度吻合，所以将理论解向上平移30mm，方便查看二者的区别。

瞬态动力学随着初始响应条件和自由伴随振动逐渐衰减，仅剩下受迫振动响应，即为谐响应计算结果，对比两者的仿真结果，在频率和幅值上，二者计算结果吻合。（虽然存在一定的误差，但随着计算时间的推移，计算结果应该会逐渐逼近）

图6. 瞬态动力学仿真解与理论解对比

图7. 瞬态动力学模拟（初始条件：位移1mm，速度0m/s；正弦载荷）

图8. 瞬态动力学模拟（初始条件：位移1mm，速度0m/s；无外部激励）

五、附录

5.1、大阻尼情形响应曲线绘图

import numpy as np
import matplotlib.pyplot as plt


# Define the constants given in the equation
k = 39.48                 # Spring constant, unit: N/m
c = 16.0                   # Damping constant, unit: N·s/m
m = 1.0                   # Quality, unit: Kg

# Intermediate variable
omega_0 = np.sqrt(k/m)                                                # Natural frequency
n = c/m*0.5


# Initial conditions
x_0 = 1.0      # The displacement at time 0, unit: m
x_dot_0 = 0.0  # The velocity at time 0, unit: m/s

# Define time array
t = np.linspace(0, 12, 1000)

# Define the function x(t)
C_1 = (x_dot_0-x_0*(-n - np.sqrt(n * n - omega_0 * omega_0))) / (2 * np.sqrt(n * n - omega_0 * omega_0))
r_1 = np.exp((-n + np.sqrt(n * n - omega_0 * omega_0)) * t)
C_2 = (-x_dot_0 + x_0*(-n + np.sqrt(n * n - omega_0 * omega_0))) / (2 * np.sqrt(n * n - omega_0 * omega_0))
r_2 = np.exp((-n - np.sqrt(n * n - omega_0 * omega_0)) * t)

x_t = C_1*r_1 + C_2*r_2


# Plotting the curve
plt.figure(figsize=(10, 6))
plt.xlim(0, 6)
plt.ylim(-1.05, 1.05)

plt.plot(t, x_t, label=r'$x(t)$', linestyle='-', linewidth=1.5)

plt.title('Free vibration response')
plt.xlabel('Time / s')
plt.ylabel('Amplitude / m')
plt.grid(True)
plt.legend()

plt.savefig('sin_curve.png', dpi=300)  # 保存图像为文件
np.savetxt('curve_data.csv', np.column_stack((t, x_t)), delimiter=',', header='x,y', comments='')  # 保存曲线数据为CSV文件
plt.show()

5.2、临界阻尼情形响应曲线绘图

import numpy as np
import matplotlib.pyplot as plt


# Define the constants given in the equation
k = 39.48                 # Spring constant, unit: N/m
c1 = 16
c2 = 12.567                # Damping constant, unit: N·s/m
m = 1.0                   # Quality, unit: Kg

# Intermediate variable
omega_0 = np.sqrt(k/m)                                                # Natural frequency
n1 = c1/m*0.5
n2 = c2/m*0.5

# Initial conditions
x_0 = 1.0      # The displacement at time 0, unit: m
x_dot_0 = 0.0  # The velocity at time 0, unit: m/s

# Define time array
t = np.linspace(0, 12, 1000)

# Define the function x(t)
C_1 = (x_dot_0-x_0*(-n1 - np.sqrt(n1 * n1 - omega_0 * omega_0))) / (2 * np.sqrt(n1 * n1 - omega_0 * omega_0))
r_1 = np.exp((-n1 + np.sqrt(n1 * n1 - omega_0 * omega_0)) * t)
C_2 = (-x_dot_0 + x_0*(-n1 + np.sqrt(n1 * n1 - omega_0 * omega_0))) / (2 * np.sqrt(n1 * n1 - omega_0 * omega_0))
r_2 = np.exp((-n1 - np.sqrt(n1 * n1 - omega_0 * omega_0)) * t)

x_t1 = C_1*r_1 + C_2*r_2
x_t2 = x_0*np.exp(-n2*t)+(x_dot_0+n2*x_0)*t*np.exp(-n2*t)


# Plotting the curve
plt.figure(figsize=(10, 6))
plt.xlim(0, 6)
plt.ylim(-1.05, 1.05)

plt.plot(t, x_t1, label=r'$x(t)$ Large damping', color='Gray', linestyle=':', linewidth=1.5)
plt.plot(t, x_t2, label=r'$x(t)$ Critical damping', linestyle='-', linewidth=1.5)


plt.title('Free vibration response')
plt.xlabel('Time / s')
plt.ylabel('Amplitude / m')
plt.grid(True)
plt.legend()

plt.savefig('sin_curve.png', dpi=300)  # 保存图像为文件
np.savetxt('curve1_data.csv', np.column_stack((t, x_t1)), delimiter=',', header='x,y', comments='')  # 保存曲线数据为CSV文件
np.savetxt('curve2_data.csv', np.column_stack((t, x_t2)), delimiter=',', header='x,y', comments='')  # 保存曲线数据为CSV文件
plt.show()

5.3、小阻尼情形响应曲线绘图

import numpy as np
import matplotlib.pyplot as plt


# Define the constants given in the equation
k = 39.48                 # Spring constant, unit: N/m
c = 2                     # Damping constant, unit: N·s/m
m = 1.0                   # Quality, unit: Kg

# Intermediate variable
omega_0 = np.sqrt(k/m)                                                # Natural frequency
n = c/m*0.5
omega_d = np.sqrt(omega_0*omega_0 - n*n)


# Initial conditions
x_0 = 1.0      # The displacement at time 0, unit: m
x_dot_0 = 0.0  # The velocity at time 0, unit: m/s

# Define time array
t = np.linspace(0, 12, 1000)

# Define the function x(t)
x_t = np.exp(-n*t)*(x_0*np.cos(omega_d*t)+((x_dot_0+n*x_0)/omega_d)*np.sin(omega_d*t))

# Plotting the curve
plt.figure(figsize=(10, 6))
plt.xlim(0, 6)
plt.ylim(-1.05, 1.05)

plt.plot(t, x_t, label=r'$x(t)$', linestyle='-', linewidth=1.5)


plt.title('Free vibration response')
plt.xlabel('Time / s')
plt.ylabel('Amplitude / m')
plt.grid(True)
plt.legend()

plt.savefig('sin_curve.png', dpi=300)  # 保存图像为文件
np.savetxt('curve1_data.csv', np.column_stack((t, x_t)), delimiter=',', header='x,y', comments='')  # 保存曲线数据为CSV文件
plt.show()

5.4、无阻尼情形响应曲线绘图

import numpy as np
import matplotlib.pyplot as plt


# Define the constants given in the equation
k = 39.48                 # Spring constant, unit: N/m
c = 0                     # Damping constant, unit: N·s/m
m = 1.0                   # Quality, unit: Kg

# Intermediate variable
omega_0 = np.sqrt(k/m)                                                # Natural frequency
n = c/m*0.5


# Initial conditions
x_0 = 1.0      # The displacement at time 0, unit: m
x_dot_0 = 0.0  # The velocity at time 0, unit: m/s

# Define time array
t = np.linspace(0, 12, 1000)

# Define the function x(t)
x_t = x_0*np.cos(omega_0*t)+(x_dot_0/omega_0)*np.sin(omega_0*t)

# Plotting the curve
plt.figure(figsize=(10, 6))
plt.xlim(0, 6)
plt.ylim(-1.05, 1.05)

plt.plot(t, x_t, label=r'$x(t)$', linestyle='-', linewidth=1.5)


plt.title('Free vibration response')
plt.xlabel('Time / s')
plt.ylabel('Amplitude / m')
plt.grid(True)
plt.legend()

plt.savefig('sin_curve.png', dpi=300)  # 保存图像为文件
np.savetxt('curve1_data.csv', np.column_stack((t, x_t)), delimiter=',', header='x,y', comments='')  # 保存曲线数据为CSV文件
plt.show()

5.5、方程各项响应曲线与全响应曲线绘图

import numpy as np
import matplotlib.pyplot as plt


# Define the constants
k = 39.48                 # Spring constant, unit: N/m
c = 2.0                   # Damping constant, unit: N·s/m
omega = 2*np.pi           # Excitation angular frequency
m = 1.0                   # Quality, unit: Kg
F_0 = 1.0                 # External force amplitude, unit: N

# Intermediate variable
zeta = c/np.sqrt(k*m)*0.5                                             # Damping ratio
omega_0 = np.sqrt(k/m)                                                # Natural frequency
omega_d = omega_0*np.sqrt(1-np.square(zeta))                          # Damped natural frequency
s = omega/omega_0
beta = 1.0/np.sqrt(np.square(1-np.square(s))+np.square(2*zeta*s))     # Amplitude amplification factor
b = ((F_0/k)*(1-s*s))/(np.square(1-s*s)+np.square(2*zeta*s))

if b > 0:
    theta = np.arctan((2*zeta*s)/(1-s*s))  # Phase angle
else:
    theta = np.arctan((2 * zeta * s) / (1 - s * s))+np.pi  # Phase angle

# Initial conditions
x_0 = 1.0      # The displacement at time 0, unit: m
x_dot_0 = 0.0  # The velocity at time 0, unit: m/s

# Define time array
t = np.linspace(0, 12, 1000)

# function x(t) = x_t1 + x_t2 + x_t3
x_t1 = np.exp(-zeta * omega_0 * t) * (x_0 * np.cos(omega_d * t) + (x_dot_0 + zeta * omega_0 * x_0) / omega_d * np.sin(omega_d * t))

x_t2 = F_0 / k * beta * np.exp(-zeta * omega_0 * t) * \
        (np.sin(theta) * np.cos(omega_d * t) + omega_0 / omega_d * (zeta * np.sin(theta) - np.cos(theta)) * np.sin(omega_d * t))

x_t3 = (F_0 / k) * beta * np.sin(omega_0 * t - theta)

x_t = x_t1 + x_t2 + x_t3




# Plotting the curve
plt.figure(figsize=(10, 6))
plt.xlim(0, 12)
plt.ylim(-0.65, 1.05)

# plt.plot(t, x_t1, label=r'$x(t)-1$', color='blue', linestyle='--', linewidth=1.5)
# plt.plot(t, x_t2, label=r'$x(t)-2$', color='green', linestyle=':', linewidth=1.5)
# plt.plot(t, x_t3, label=r'$x(t)-3$', color='gray', linestyle='-', linewidth=1.5)
plt.plot(t, x_t, label=r'$x(t)$', linestyle='-', linewidth=1.5)

# plt.title('Response curves of different parts')
plt.title('Full response curve')

plt.xlabel('Time / s')
plt.ylabel('Amplitude / m')
plt.grid(True)
plt.legend()

plt.savefig('sin_curve.png', dpi=300)  # 保存图像为文件
np.savetxt('curve_data.csv', np.column_stack((t, x_t)), delimiter=',', header='x,y', comments='')  # 保存曲线数据为CSV文件
plt.show()