基于频率响应函数的动力学模型修正方法研究

摘要

概述了国内外动力学模型修正技术的研究状况，研究了近些年发展起来的基于频响函数的动力学模型修正方法；利用航天器振动试验测量所得的频响函数，从理论上介绍了频响函数残差法、设计参数型频响函数法和摄动型频响函数法三种基于频响函数的动力学模型修正方法，为动力学模型修正技术的发展提供参考。

正文

《1 前言》

1 前言

随着航天器型号任务的增多, 为了降低成本, 缩短研制周期, 要求对航天器的动力学特性进行预示, 有限元模型修正技术因此变得非常重要。

近年来有限元法已发展成为结构分析的重要手段, 已经开发了行之有效的大型结构分析软件, 如ANSYS, NASTRAN等。有限元分析速度快, 设计周期短, 与结构动力学试验相比, 费用较低;使用有限元模型进行结构设计, 可以计算各种载荷、各种边界条件下的动力响应。但是, 经验表明, 拥有好的分析软件并不能保证算出正确的结果;在有好的分析软件前提下, 分析的正确性主要取决于有限元模型的好坏。而有限元建模与实际对象比较往往存在误差, 如边界条件和连接条件的简化、几何模型和本构关系的不准确, 系统阻尼必须人为引入等^[1]。对于复杂的实际结构, 特别是航天器结构, 有限元模型预示与试验结果之间往往存在明显误差, 为了消除这种误差, 有必要进行有限元模型修正。其次, 现代航天器设计要求平台化、结构化、模块化, 尽量缩短研制周期, 以适应空间市场的快速发展。在进一步进行适当试验的基础上, 在利用现有平台进行改进设计的情况下, 采用模型修正技术有可能省掉一些大型结构试验, 根据需要做适当数量的零部件级试验进行验证和修正, 以使修正后的有限元模型代表真实的结构模型。因此, 结构动力学模型修正技术在现代航天器结构设计中具有重要意义。

对航天器进行振动试验得到的数据往往是航天器不同测点的响应数据, 为了使航天器动力学模型修正技术具有更实际的意义, 必须研究基于振动响应或频响函数的动力学模型修正技术, 即利用测量的振动响应数据来调整有限元模型的参数, 使调整后的有限元模型的频响特性更为精确。据此, 笔者利用航天器力学振动试验和有限元计算得到的频响函数, 详细介绍了频响函数残差法、设计参数型频响函数法和摄动型频响函数法三种基于频响函数的动力学模型修正方法, 为航天器动力学模型修正技术的发展提供有益的参考。

《2 动力学模型修正方法研究状况》

2 动力学模型修正方法研究状况

20世纪70年代以来, 国内外研究人员一直在对试验/分析相关性和有限元模型修正技术进行研究, 并提出了许多相关性分析方法和模型修正方法。早期, 人们将系统响应的试验/分析相关性当作结构系统辨识方法来研究^[2,3], 这些方法通过调整数学模型的参数来改进分析模型, 如Collins等首先提出了一种使用测量频率和振型对分析模型进行刚度和质量特性修正的统计参数估计方法^[2], 这种方法通过修正物理结构参数和设计变量, 间接地改进分析模型的刚度和质量矩阵。后来, Hasselman和Chrostowski在大量试验的基础上使用一种统计参数估计方法来修正物理结构的参数^[4], 并编制了计算机程序代码SSID。其他模型修正统计参数估计方法的程序代码还有Dascotte的SYSTUNE^[5]。Chen和Garba提出了一种利用模态试验得到的模态数据矩阵和特征值矩阵来修改结构的物理参数的矩阵摄动法, 并对某大型空间结构进行了在轨损伤评价^[6,7,8]。Muralidhar和Thomas等采用试验测量和分析所得的模态数据对有限元模型进行修正^[9]。Berman和Flannelly提出了一种使用不完全测量模态对一个线性的、离散的结构模型进行参数识别的模型修正分析方法^[10]。随后, Baruch和Bar Itzhack提出了一种用于获得测量振型的正交性模态的方法^[11], 用于改进分析模型的刚度矩阵。Berman和Nagy应用约束最小化理论, 采用测量振型和频率来改进结构分析模型的质量和刚度矩阵^[12], 该方法假定测量振型是准确的, 并应用标准拉格朗日乘子方法解决优化问题。Kabe应用结构相关性信息来调整不正确的分析刚度矩阵^[13]。Wei导出一种元素修正方法来修改动力学模型^[14,15], 应用结构相关性原理同时修正结构刚度和质量矩阵。Chou, O´Callahan和Wu应用在模型空间中结构单元的相关性原理编制了一个用于确定单元质量和刚度矩阵误差位置的程序^[16]。Srinivas与Kao等采用动力参数灵敏度对某一行星际飞船进行分析/试验相关性分析^[17], 用自编的误差定位程序和参数估计方法程序与NASTRAN软件相结合, 对飞船有限元模型进行误差定位和参数估计。

国内研究人员在对动力学模型修正技术研究的初期, 质量矩阵和刚度矩阵的修正都是基于正交性条件求解的。考虑到质量矩阵和刚度矩阵的相关性, 张德文发展的方法是质量矩阵用正交性条件修正和刚度矩阵用特征方程修正^[18], 从算法上, 张德文还定义了一种广义逆, 用正交性条件和特征方程进行模型修正^[19]。此外, 还有向锦武等提出了基于矩阵型模型修正的两步法^[20]:先通过误差矩阵判断模型误差的范围, 然后在该范围内确定修正量。李书等提出了一种带约束条件的广义特征值模型修正方法^[21], 仅需少量的试验数据, 其有效性通过一个简单结构的数值计算得到验证。张令弥等提出了一种设计灵敏度分析的迭代模态法^[22], 可以提高特征向量对设计参数导数的计算效率, 从理论上证明了经典模态法和修改模态法只是迭代模态法的特殊情况。张景绘系统地介绍了一些动力学模型修正的方法^[23]。

模型修正方法可以划分为模态法和频响函数法。从试验特征值和特征向量来修正有限元模型参数的方法称模态法, 其关键在于试验模态参数识别的精度;基于频响函数的模型修正是近些年发展起来的, 如徐张明等利用试验测试和有限元模型计算得到的频响函数^[24,25], 推导出了一种基于频响函数的灵敏度分析模型修正方法。秦仙蓉等详细、深入地研究了基于灵敏度分析的模型修正方法^[26], 并从理论上讨论了一种频响函数残差法。Hemez等分析了应用频响函数进行动力学模型修正的难点^[27]。Ting提出了一种改进的计算频响函数灵敏度的方法^[28], 用于研究频响函数模型修正技术等。频响函数法回避了试验模态参数识别这个步骤, 直接利用有限元计算和测量得到的频响函数进行模型修正, 其优点是:克服了模态法需要振型一一对应的缺点, 并且修正计算的频率范围很宽, 采用适当数目频率点上的频响函数, 便可求解问题^[24]。

《3 基于频响函数的模型修正方法》

3 基于频响函数的模型修正方法

《3.1 频响函数残差法[26]》

3.1 频响函数残差法[26]

n个自由度线性时不变系统的运动可以由以下微分方程描述:

$Μ \ddot{x} + C \dot{x} + Κ x = f (1)$ $Μ \ddot{x} + C \dot{x} + Κ x = f (1)$

式中K, M, C, 依次为系统的刚度矩阵、质量矩阵和阻尼矩阵;x, f是位移响应向量与输入力激励向量。方程式 (1) 两边同时进行傅里叶变换可得:

$\begin{array}{l} Ζ (ω) X (ω) = F (ω) ‚ \\ Ζ (ω) = - ω^{2} Μ + j ω C + Κ (2) \end{array}$ $\begin{array}{l} Ζ (ω) X (ω) = F (ω) ‚ \\ Ζ (ω) = - ω^{2} Μ + j ω C + Κ (2) \end{array}$

式中j, ω分别代表复数单位和频率, X (ω) , F (ω) 是位移响应与输入力激励的傅里叶变换, Z (ω) 为结构位移阻抗矩阵 (动刚度矩阵) , 其逆矩阵即结构的频响函数矩阵H (ω) :

《图1》

则结构位移响应可以根据式 (2) 表述为

$X (ω) = Η (ω) F (ω) (4)$ $X (ω) = Η (ω) F (ω) (4)$

若输入力激励为一个作用于结构第j个自由度的单位力, 则上式可以简化为

$X = {Η}_{j} (5)$ $X = {Η}_{j} (5)$

这里{H}_j代表频响矩阵的第j列。

定义残差项为结构位移响应的实验测试值与理论预测值的差值, 有:

$R (p) = X_{t} - X_{A} (p) (6)$ $R (p) = X_{t} - X_{A} (p) (6)$

这里p代表设计参数, X_t, X_A表示位移响应的试验值和理论值。将式 (4) 、式 (5) 代入式 (6) 可得

$\begin{array}{l} R (p) = {Η_{t}}_{j} - Η_{A} (p) {Ι}_{j} = \\ {Η_{t}}_{j} - Η_{A} (p) ； (7) \end{array}$ $\begin{array}{l} R (p) = {Η_{t}}_{j} - Η_{A} (p) {Ι}_{j} = \\ {Η_{t}}_{j} - Η_{A} (p) ； (7) \end{array}$

实测频响函数一般为复数, 所以在这种残差定义下, 修正问题的目标函数也是一个复数。在一系列选定的频率点, 如果结构位移响应的理论预测值与试验测试值具有良好的相关性, 则残差项R (p) 的各个分量应该小于用户给定的数值, 或者残差项的模‖R (p) ‖在这一系列频率点位置得到最小化。这就是基于结构频率响应函数的模型修正技术的基本思路。

以上定义的残差是结构位移响应残差。这种表达虽然很直观, 但由于频响函数矩阵的预测值H_A (ω) 在给定频率点相对于设计参数p的变化不光滑, 并且对低阻尼系统H_A (ω) 在共振点附近也不连续, 实际过程一般不直接利用位移响应残差。

为此, Link提出给位移响应形式的残差左乘阻抗矩阵^[29], 利用结构的力响应信息来修正模型。定义残差项为输入力激励的试验测试值与理论预测值之间的差异:

$R (p) = F_{E} - F_{A} (ω, p) (8)$ $R (p) = F_{E} - F_{A} (ω, p) (8)$

由于F_A (p) =Z_A (ω, p) X_A (p) , 式 (8) 可以改写为

$R (p) = F_{E} - Ζ_{A} (ω, p) X_{A} (p) (9)$ $R (p) = F_{E} - Ζ_{A} (ω, p) X_{A} (p) (9)$

与上面类似, 在结构第j个自由度施加一个单位力, 根据式 (5) 有

$R (p) = {Ι}_{j} - Ζ_{A} (ω, p) Η_{E}_{j} (10)$ $R (p) = {Ι}_{j} - Ζ_{A} (ω, p) Η_{E}_{j} (10)$

式中{I}_j代表第j个分量为1, 其余分量为0的列向量。式 (10) 代表一系列非线性方程, 可以利用一阶泰勒公式线性化。

设结构动刚度矩阵在p₀的一阶泰勒展开式为

$Ζ_{A} (ω, p) = Ζ_{A}^{0} + \sum_{i} \frac{\partial Ζ_{A}}{\partial p_{i}} Δ p_{i} (11)$ $Ζ_{A} (ω, p) = Ζ_{A}^{0} + \sum_{i} \frac{\partial Ζ_{A}}{\partial p_{i}} Δ p_{i} (11)$

式中Z $_{A}^{0}$ $_{A}^{0}$ =Z_A (ω, p₀) , Δp=[Δp₁Δp₂ … Δp_{N_p}]^T是设计参数变化量, N_p为设计参数个数。这样式 (10) 可以最终表示为

$R (p) = {Ι}_{j} - Ζ_{A}^{0} {Η_{E}}_{j} - (\sum_{i} \frac{\partial Ζ_{A}}{\partial p_{i}} Δ p_{i}) {Η_{E}}_{j} (12)$ $R (p) = {Ι}_{j} - Ζ_{A}^{0} {Η_{E}}_{j} - (\sum_{i} \frac{\partial Ζ_{A}}{\partial p_{i}} Δ p_{i}) {Η_{E}}_{j} (12)$

则模型修正问题可以通过求解下面的最小化问题得到解决:

$\begin{array}{l} \min_{p} ∥ R (p) ∥_{2}^{2} (13) \\ s . t V L B \leq p \leq V U B \end{array}$ $\begin{array}{l} \min_{p} ∥ R (p) ∥_{2}^{2} (13) \\ s . t V L B \leq p \leq V U B \end{array}$

式中VLB, VUB是设计参数p变化的上下限。

《3.2 设计参数型频响函数法》

3.2 设计参数型频响函数法

该法将实际系统的参数矩阵表示为计算模型的参数矩阵与摄动量之和, 而摄动量为参数的一阶摄动, 在频率域内利用试验测量加速度响应来求设计参数的修正量^[23]。方法如下:

考虑n自由度粘性阻尼系统, 在简谐激励下, 频域内输入和输出的关系为

$X (j ω) = Η (j ω) F (j ω) (14)$ $X (j ω) = Η (j ω) F (j ω) (14)$

式中X (jω) 为稳态加速度响应;H (jω) 为加速度频响函数;F (jω) 为简谐激励。

设有限元计算模型和试验模型的加速度频率响应函数矩阵为

$\begin{array}{l} \begin{array}{l} Η_{A} (j ω) = - \frac{ω^{2}}{(- ω^{2} Μ_{A} + j ω C_{A} + Κ_{A})} = \\ - (- Μ_{A} + j C_{A} / ω + Κ_{A} / ω^{2})^{- 1} (15) \end{array} \\ \begin{array}{l} Η_{t} (j ω) = - \frac{ω^{2}}{(- ω^{2} Μ_{t} + j ω C_{t} + Κ_{t})} = \\ - (- Μ_{t} + j C_{t} / ω + Κ_{t} / ω^{2})^{- 1} (16) \end{array} \end{array}$ $\begin{array}{l} \begin{array}{l} Η_{A} (j ω) = - \frac{ω^{2}}{(- ω^{2} Μ_{A} + j ω C_{A} + Κ_{A})} = \\ - (- Μ_{A} + j C_{A} / ω + Κ_{A} / ω^{2})^{- 1} (15) \end{array} \\ \begin{array}{l} Η_{t} (j ω) = - \frac{ω^{2}}{(- ω^{2} Μ_{t} + j ω C_{t} + Κ_{t})} = \\ - (- Μ_{t} + j C_{t} / ω + Κ_{t} / ω^{2})^{- 1} (16) \end{array} \end{array}$

令计算模型的初始矩阵为M_A, C_A, K_A, 试验模型的初始矩阵为M_t, C_t, K_t, 对计算模型的M_A, C_A, K_A进行修正, 引入一个小摄动量, 有

$\begin{array}{l} Μ_{t} = Μ_{A} + Δ Μ ‚ \\ C_{t} = C_{A} + Δ C ‚ \\ Κ_{t} = Κ_{A} + Δ Κ (17) \end{array}$ $\begin{array}{l} Μ_{t} = Μ_{A} + Δ Μ ‚ \\ C_{t} = C_{A} + Δ C ‚ \\ Κ_{t} = Κ_{A} + Δ Κ (17) \end{array}$

假设a_i为某个设计参数, Δa_i为修正量, 即变化量ΔM, ΔC, ΔK近似地表示为修正参数的一阶摄动, 有

《图2》

设X_A (jω) 为计算模型在力F_A (jω) 作用下的响应, 将式 (15) 代入式 (14) 得

$- (- Μ_{A} + j C_{A} / ω + Κ_{A} / ω^{2}) X_{A} (j ω) = F_{A} (j ω) (19)$ $- (- Μ_{A} + j C_{A} / ω + Κ_{A} / ω^{2}) X_{A} (j ω) = F_{A} (j ω) (19)$

试验测量响应Xt (jω) 在力F_t (jω) 作用下的响应为

$- (- Μ_{t} + j C_{t} / ω + Κ_{t} / ω^{2}) X_{t} (j ω) = F_{t} (j ω) (20)$ $- (- Μ_{t} + j C_{t} / ω + Κ_{t} / ω^{2}) X_{t} (j ω) = F_{t} (j ω) (20)$

由于计算模型的作用力与试验测量的作用力相同, 即

$F_{A} (j ω) = F_{t} (j ω) (21)$ $F_{A} (j ω) = F_{t} (j ω) (21)$

将式 (17) 、式 (19) 、式 (20) 代入式 (21) 可得

$\begin{array}{l} (- Μ_{A} + j C_{A} / ω + Κ_{A} / ω^{2}) (X_{A} - X_{t}) = \\ (- Δ Μ + j Δ C / ω + Δ Κ / ω^{2}) X_{t} (22) \end{array}$ $\begin{array}{l} (- Μ_{A} + j C_{A} / ω + Κ_{A} / ω^{2}) (X_{A} - X_{t}) = \\ (- Δ Μ + j Δ C / ω + Δ Κ / ω^{2}) X_{t} (22) \end{array}$

由于

《图3》

有

$\begin{array}{l} Η_{A} (j ω) (- Δ Μ + j Δ C / ω + Δ Κ / ω^{2}) X_{t} = \\ (X_{t} - X_{A}) (23) \end{array}$ $\begin{array}{l} Η_{A} (j ω) (- Δ Μ + j Δ C / ω + Δ Κ / ω^{2}) X_{t} = \\ (X_{t} - X_{A}) (23) \end{array}$

再将式 (18) 代入式 (23) 中, 得

$\begin{array}{l} \sum_{i} [Η_{A} (j ω) (- \frac{\partial Μ}{\partial a_{i}} + j \frac{\partial C}{\partial a_{i}} \frac{1}{ω_{A}} + \frac{\partial Κ}{\partial a_{i}} \frac{1}{ω_{A}^{2}}) \cdot \\ X_{t} (j ω)] Δ a_{i} = X_{t} (j ω) - X_{A} (j ω) (24) \end{array}$ $\begin{array}{l} \sum_{i} [Η_{A} (j ω) (- \frac{\partial Μ}{\partial a_{i}} + j \frac{\partial C}{\partial a_{i}} \frac{1}{ω_{A}} + \frac{\partial Κ}{\partial a_{i}} \frac{1}{ω_{A}^{2}}) \cdot \\ X_{t} (j ω)] Δ a_{i} = X_{t} (j ω) - X_{A} (j ω) (24) \end{array}$

由此得到一个复数形式的关于未知修正参数Δa_i的线性方程。

《3.3 摄动型频响函数法》

3.3 摄动型频响函数法

该法通过引入一个摄动量, 将质量、阻尼、刚度的修正量局部化, 并利用位置矩阵进行一系列的矩阵运算, 建立被修正参数的估计表达式, 由曲线拟合得出质量、阻尼、刚度的局部修正量^[23]。

将式 (17) 代入式 (16) 可得

$\begin{array}{l} Η_{t} (j ω) = - (- Μ_{A} + j C_{A} / ω + Κ_{A} / ω^{2} - \\ Δ Μ + j Δ C / ω + Δ Κ / ω^{2})^{- 1} = (Ζ_{A} + D)^{- 1} \end{array} (25)$ $\begin{array}{l} Η_{t} (j ω) = - (- Μ_{A} + j C_{A} / ω + Κ_{A} / ω^{2} - \\ Δ Μ + j Δ C / ω + Δ Κ / ω^{2})^{- 1} = (Ζ_{A} + D)^{- 1} \end{array} (25)$

式中Z_A=- (-M_A+jC_A/ω+K_A/ω²) =

H $_{A}^{- 1}$ $_{A}^{- 1}$ (jω) (26)

D=- (-ΔM+jΔC/ω+ΔK/ω²) (27)

当结构比较复杂时, 只能测量得到部分信息。在这种不完备信息的情况下, ΔM, ΔC, ΔK也没有唯一解。这里将质量、阻尼、刚度的修正量局部化, 利用位置矩阵将修正量集中体现在测量自由度上。将式 (17) 表示为

$\begin{array}{l} Δ Μ = Ρ m Ρ^{Τ} ‚ \\ Δ [C X 2] C = Ρ c Ρ^{Τ} ‚ \\ Δ [C X 2] Κ = Ρ k Ρ^{Τ} (28) \end{array}$ $\begin{array}{l} Δ Μ = Ρ m Ρ^{Τ} ‚ \\ Δ [C X 2] C = Ρ c Ρ^{Τ} ‚ \\ Δ [C X 2] Κ = Ρ k Ρ^{Τ} (28) \end{array}$

式中m, c, k均为r×r阶的实对称矩阵, P∈R^n×r称之为位置矩阵 (r为测量自由度数目) , 它是由n阶单位阵划去非测量自由度对应列形成的矩阵。

将式 (28) 代入式 (27) 得

$\begin{array}{l} D = - (- Ρ m Ρ^{Τ} + j Ρ c Ρ^{Τ} / ω + Ρ k Ρ^{Τ} / ω^{2}) = \\ Ρ \tilde{D} Ρ^{Τ} (29) \end{array}$ $\begin{array}{l} D = - (- Ρ m Ρ^{Τ} + j Ρ c Ρ^{Τ} / ω + Ρ k Ρ^{Τ} / ω^{2}) = \\ Ρ \tilde{D} Ρ^{Τ} (29) \end{array}$

$式中 \tilde{D} = - (- m + j c / ω + k / ω^{2}) (30)$ $式中 \tilde{D} = - (- m + j c / ω + k / ω^{2}) (30)$

将式 (29) 代入式 (25) 得

$Ζ_{A} Η_{t} + Ρ \tilde{D} Ρ^{Τ} Η_{t} = Ι (31)$ $Ζ_{A} Η_{t} + Ρ \tilde{D} Ρ^{Τ} Η_{t} = Ι (31)$

将式 (31) 左乘以P^TH_A, 右乘以P, 得

$Ρ^{Τ} Η_{A} Ζ_{A} Η_{t} Ρ + Ρ^{Τ} Η_{A} Ρ \tilde{D} Ρ^{Τ} Η_{t} Ρ = Ρ^{Τ} Η_{A} Ρ (32)$ $Ρ^{Τ} Η_{A} Ζ_{A} Η_{t} Ρ + Ρ^{Τ} Η_{A} Ρ \tilde{D} Ρ^{Τ} Η_{t} Ρ = Ρ^{Τ} Η_{A} Ρ (32)$

将式 (26) 代入式 (32) , 有

$Ρ^{Τ} Η_{t} Ρ + Ρ^{Τ} Η_{A} Ρ \tilde{D} Ρ^{Τ} Η_{t} Ρ = Ρ^{Τ} Η_{A} Ρ (33)$ $Ρ^{Τ} Η_{t} Ρ + Ρ^{Τ} Η_{A} Ρ \tilde{D} Ρ^{Τ} Η_{t} Ρ = Ρ^{Τ} Η_{A} Ρ (33)$

令

$\begin{array}{l} Η_{f} = Ρ^{Τ} Η_{A} Ρ (34) \\ Η_{m} = Ρ^{Τ} Η_{t} Ρ (35) \end{array}$ $\begin{array}{l} Η_{f} = Ρ^{Τ} Η_{A} Ρ (34) \\ Η_{m} = Ρ^{Τ} Η_{t} Ρ (35) \end{array}$

代入式 (33) 得

$Η_{m} + Η_{f} \tilde{D} Η_{m} = Η_{f} (36)$ $Η_{m} + Η_{f} \tilde{D} Η_{m} = Η_{f} (36)$

所以

《图4》

由式 (30) 有

《图5》

式中m_ij, c_ij, k_ij分别为m, c, k的元素, $\tilde{D}$ $\tilde{D}$ _ij为 $\tilde{D}$ $\tilde{D}$ 的元素。根据结构振动系统要求的频率范围, 利用曲线拟合法, 由式 (38) 可以估计得到每个m_ij, c_ij, k_ij的值, 然后由式 (28) 和式 (17) 得到修正后的系统矩阵。

《4 结论》

4 结论

尽管动力学模型修正技术还在发展, 并有许多理论方法问题尚需研究, 但是工程实践表明, 模型修正技术在航天器结构设计中具有重要意义。它不但可以提高航天器的动力学特性预示精度, 利用先验试验数据有时可以省略结构星的研制和试验, 大量节约研制经费, 缩短研制周期和提高设计质量。

动力学模型修正技术是一门有着广阔前景的新技术, 基于频响函数的模型修正方法能够直接利用力学振动试验测量所得的频响函数矩阵, 具有较大的工程实用价值而备受人们的青睐。目前频响函数修正方法还处在探索和研究阶段, 问题较多, 如振动测试数据与程序接口、修正频段的选择、参加修正的频率点数的确定。在对振动试验测量所得的频响函数的预处理中, 很难判断哪些是正常的, 哪些是异常的。实现模型修正软件开发和解决编程语言如C++与有限元分析软件如NASTRAN软件的接口难度较大。所以目前的频响函数修正方法在大型复杂结构上应用, 尤其是航天器结构设计上应用, 尚有一定的距离, 许多问题有待研究解决。

展示更多