CN106384384A - 一种三维产品模型的形状优化算法 - Google Patents

Abstract

本发明提供了一种三维产品模型的形状优化算法，包括以下步骤：首先利用反求方法或者尺寸测量方法，构建产品模型包含点、线、面特征的特征框架模型，添加特征元素之间的约束，构成参数化的特征框架，通过多种驱动方式且维护特征之间的约束以实现特征框架模型的修改。再针对该特征框架，利用体参数化的理论得到可用于等几何分析的体参数化模型。同时将该特征框架的特征尺寸作为产品尺寸或者形状优化的参数，推导产品物理性能与这些设计参数的敏度矩阵，通过优化这些参数以实现对产品模型的优化设计。

Description

一种三维产品模型的形状优化算法

技术领域

本发明属于三维产品模型形状优化设计领域，主要涉及基于产品模型特征框架的形状优化算法。

背景技术

在传统形状优化中，通常将有限单元的结点作为设计变量。但是在这种方法中，由于各结点独立的变动，导致优化边界常常是不规则的，设计变量过多导致优化效率降低，设计模型与分析模型之间在优化迭代中需要不断进行数据转换，网格划分费时费力。而这些缺点都是由于设计模型与分析模型采用不同的数学语言引起对的，即几何建模中采用样条基函数，而在分析模型中采用拉格朗日基函数或埃尔米特多项式。

发明内容

本发明是为解决上述问题而进行的，通过建立参数化的特征框架模型作为产品的优化对象，将特征框架模型的尺寸作为设计变量。利用体插值理论通过特征框架模型生成参数化的控制点坐标，推导目标函数及约束对设计变量的灵敏度矩阵进行特征框架模型的尺寸优化。由于等几何分析所采用的B样条基函数具有高阶连续性，其导数具有准确的数学表达式，可以得到灵敏度的全解析式，进而可以高效准确地进行优化设计。

本发明所提供的三维产品模型的形状优化算法，包括以下四个步骤：

(一)建立特征框架模型

利用模型扫描得到点云数据，提取点、线、面特征信息，或者利用模型测量和尺寸设计，通过这几种方法得到自由形体模型；根据设计对象的特征和设计意图，选取部分特征参数作为优化设计参数，根据这些参数和其他一些已知参数，构建出产品模型的特征曲线；给所述特征曲线添加约束条件得到参数化的特征框架，通过特征框架尺寸驱动产品模型的变形。

(二)建立分析模型

采用一定方法获取特征框架模型的表面控制点和实体内部控制点，并且建立优化设计参数和控制点之间的关联。

(三)确定优化算法

给特征框架模型施加边界条件或约束条件，利用等几何分析的方法得到每个控制点的物性参数；利用控制点和优化设计参数的关系，从而构建出优化设计参数和物性参数的关系；根据该关系求取物性参数针对于优化设计参数的导数关系，即得到灵敏度矩阵；对灵敏度矩阵进行求解，确定优化算法。

(四)利用选用的优化算法，依次对每个优化设计参数和对应的特征框架模型进行更新，直到满足给定的迭代终止条件。

发明作用与效果

本发明提供的基于特征框架的优化方法与以控制点坐标为设计变量或者以有限元网格顶点为优化变量的传统优化方法相比，具有明显优势：在设计理念上，设计变量更贴近设计者的意图，符合人们的设计需要，重点保证和优化关键部位尺寸，利用特征框架模型信息实现对模型总体形状和细节形状的控制，使得用户能通过少数特征尺寸参数编辑产品模型，避免传统形状优化以有限元节点为设计变量带来的复杂的造型操作，将参数化建模与优化真正统一起来；在计算效率上，因为大大减少了设计变量，带来计算效率的提高；在设计流程上，设计模型、分析模型、优化模型之间的交互更加便捷，大大降低产品开发的周期。

附图说明

图1是本发明中的三维产品模型的形状优化算法的流程图；

图2是本发明所采用的离散孔斯插值法的示意图；

图3是本发明实验例一中的带脊悬臂梁的特征框架模型的示意图；

图4是本发明实验例一中的带脊悬臂梁的优化历史图；

图5是本发明实验例二中的疏浚绞刀臂的特征框架模型的示意图；

图6是本发明实验例二中的疏浚绞刀臂的优化历史图。

具体实施方式

以下结合附图来说明本发明的具体实施方式。

图1是本发明实施例中的三维产品模型的形状优化算法的流程图。

如图1所示，本发明中的三维产品模型的形状优化的一般算法包括以下部分：

步骤1，构建特征框架模型

步骤1-1，利用模型扫描得到点云数据，提取点、线、面特征信息，或者利用模型测量和尺寸设计，得到自由形体模型的特征框架；

步骤1-2，根据设计对象的特征和设计意图，选取部分特征参数作为优化设计参数，根据这些参数和其他一些已知参数，构建出产品模型的特征曲线；

步骤1-3，给所述特征曲线添加约束条件得到参数化的特征框架，通过特征框架尺寸驱动产品模型的变形，得到特征框架模型；

步骤2，建立分析模型

从特征框架模型到分析模型的生成，采用经典的孔斯插值方法，可以分为两步。第一步为通过特征曲线生成边界面所有控制点，第二步为通过所有表面控制点生成实体内部所有控制点。下面详细描述插值方法：

步骤2-1，如图2(a)所示，已知四条曲面边界线上的控制点，通过孔斯曲面插值，得到曲面上所有控制点坐标。设某边界面控制点为P(i,j)，n为U方向控制点个数，m为V方向控制点个数。其内部控制点的下标范围为1<i<m，1<j<n。

令u₀＝(i-1)/(m-1)，u₁＝1-u0；v₀＝(j-1)/(n-1)，v₁＝1-v0，其中1<i<m，1<j<n。则边界面控制点P(i,j)通过式(1)得到：

$P_{i,j}=u1*p_{0,j}+u0*p_{m,j}+v1*p_{i,0}+v0*p_{\underset{-}{i,n}}-\&lsqb;\begin{array}{cc}u1& u0\end{array}\&rsqb;*\left[\begin{array}{cc}{p}_{\underset{-}{i,n}}& p_{\underset{-}{i,n}}\\ p_{\underset{-}{i,n}}& p_{\underset{-}{i,n}}\end{array}\right]\left[\begin{array}{c}v1\\ v0\end{array}\right]---\left(1\right)$

步骤2-2，如图2(b)所示，已知实体边界面上的控制点，通过孔斯体插值，得到实体内部所有控制点坐标。设某边界面控制点为P(i,j,k)，l为U方向控制点数目，m为V方向控制点数目，n为W方向上的控制点数目。其内部控制点的下标范围为1<i<l,1<j<m,1<k<n。

令u₀＝(i-1)/(m-1)，u₁＝1-u₀；v₀＝(j-1)/(n-1)，v₁＝1-v₀；w₀＝(k-1)/(n-1)，w₁＝1-w₀，其中1<i<l,1<j<m,1<k<n。则实体内部控制点P(i,j,k)通过式

(2)得到：

$\begin{array}{c}{P}_{i,j,k}=u1*p_{0,j,k}+u0*p_{l,j,k}+v1*p_{i,0,k}+v0*p_{\underset{-}{i,m,k}}+w1*p_{i,j,0}+w0*p_{\underset{-}{i,j,n}}-\\ \left[\begin{array}{cc}u1& u0\end{array}\right]\left[\begin{array}{cc}{p}_{\underset{-}{0,0,k}}& p_{\underset{-}{0,m,k}}\\ p_{\underset{-}{l,0,k}}& p_{\underset{-}{l,m,k}}\end{array}\right]\left[\begin{array}{c}v1\\ v0\end{array}\right]-\left[\begin{array}{cc}v1& v0\end{array}\right]\left[\begin{array}{cc}pi{,}_{\underset{-}{0,0}}& p_{\underset{-}{i,0,n}}\\ p_{\underset{-}{i,m,0}}& p_{\underset{-}{i,m,n}}\end{array}\right]\left[\begin{array}{c}w1\\ w0\end{array}\right]-\left[\begin{array}{cc}w1& w0\end{array}\right]\left[\begin{array}{cc}{p}_{\underset{-}{0,j,0}}& p_{\underset{-}{l,j,0}}\\ p_{\underset{-}{0,j,n}}& p_{\underset{-}{l,j,n}}\end{array}\right]\left[\begin{array}{c}u1\\ u0\end{array}\right]+\\ w1*\&lsqb;\left[\begin{array}{cc}u1& u0\end{array}\right]\left[\begin{array}{cc}{p}_{\underset{-}{0,0,0}}& p_{\underset{-}{0,m,0}}\\ p_{\underset{-}{l,0,0}}& p_{\underset{-}{l,m,0}}\end{array}\right]\left[\begin{array}{c}v1\\ v0\end{array}\right]\&rsqb;+w0*\&lsqb;\left[\begin{array}{cc}u1& u0\end{array}\right]\left[\begin{array}{cc}{p}_{\underset{-}{0,0,n}}& p_{\underset{-}{0,m,n}}\\ p_{\underset{-}{l,0,n}}& p_{\underset{-}{l,m,n}}\end{array}\right]\left[\begin{array}{c}v1\\ v0\end{array}\right]\&rsqb;\end{array}---\left(2\right)$

步骤3，建立控制点坐标对优化设计参数的灵敏度矩阵

步骤3-1，向特征框架模型施加边界条件或约束条件，利用等几何分析的方法得到每个控制点的位移和受力等物性参数；

步骤3-2，利用控制点和优化设计参数的关系，从而构建出优化设计参数和物性参数的关系；

步骤3-3，根据该关系求取物性参数针对于优化设计参数的导数关系，即得到灵敏度矩阵，并对其进行求解，得到优化算法：

首先，求解目标方程和约束方程对设计变量灵敏度，在本发明中，优化目标为模型刚度最大，约束为模型体积小于特定值。

下面求解目标方程和约束方程对设计变量(优化设计参数)的导数。

(1)三维产品模型的结构柔度对优化设计参数的导数为：

$\frac{d\left(f^{T}u\right)}{{\mathrm{d\&alpha;}}_i}=f^T\frac{\&part;u}{\&part;{\mathrm{\&alpha;}}_i}=-f^T{K}^{-1}\frac{\&part;K}{\&part;{\mathrm{\&alpha;}}_i}u=-u^T\frac{\&part;K}{\&part;{\mathrm{\&alpha;}}_i}u---\left(3\right)$

(2)三维产品模型的体积或面积对优化设计参数的导数为：

$\frac{dS}{{\mathrm{d\&alpha;}}_i}=\underset{e=1}{\overset{n_{el}}{\&Sigma;}}\frac{d\left|J\right|}{{\mathrm{d\&alpha;}}_i}d\widetilde{\mathrm{\&Omega;}}---\left(4\right)$

式中，f为载荷列阵；u为位移列阵；α_i为优化设计参数；K为三维产品模型的整体刚度矩阵；S为三维产品模型的体积或面积；|J|为雅克比行列式；Ω为分析域。

其次，对灵敏度矩阵进行求解，

灵敏度分析首先要求离散控制方程式(13)对设计变量的导数。而整体刚度矩阵K是由单元刚度K_e矩阵组装而来，因此转化为求单元刚度矩阵K_e对设计变量α_i的导数问题：

$\frac{\&part;K_e}{\&part;{\mathrm{\&alpha;}}_i}=\underset{{\widetilde{\mathrm{\&Omega;}}}_e}{\&Integral;}\frac{\&part;B^T}{\&part;\mathrm{\&alpha;}}DB\left|J\right|d{\widetilde{\mathrm{\&Omega;}}}_e+\underset{{\widetilde{\mathrm{\&Omega;}}}_e}{\&Integral;}{B}^{T}D\frac{\&part;B}{\&part;\mathrm{\&alpha;}}\left|J\right|d{\widetilde{\mathrm{\&Omega;}}}_{ee}+\underset{{\widetilde{\mathrm{\&Omega;}}}_e}{\&Integral;}{B}^{T}DB\frac{\&part;\left|J\right|}{\&part;\mathrm{\&alpha;}}d{\widetilde{\mathrm{\&Omega;}}}_e---\left(5\right)$

根据式(3)和(4)，需要计算B及|J|对设计变量x的导数，首先计算B对设计变量x(和设计变量α_i等同)的导数：

B＝[B₁ B₂ ... B_i ... B_n] (6)

$B_i=\left[\begin{array}{ccc}\frac{\&part;N_i}{\&part;x}& 0& 0\\ 0& \frac{\&part;N_i}{\&part;y}& 0\\ 0& 0& \frac{\&part;N_i}{\&part;z}\\ \frac{\&part;N_i}{\&part;y}& \frac{\&part;N_i}{\&part;x}& 0\\ 0& \frac{\&part;N_i}{\&part;z}& \frac{\&part;N_i}{\&part;y}\\ \frac{\&part;N_i}{\&part;z}& 0& \frac{\&part;N_i}{\&part;x}\end{array}\right]---\left(7\right)$

定义两个矩阵

$T=\left[\begin{array}{cccc}\frac{\&part;N_1}{\&part;x}& \frac{\&part;N_2}{\&part;x}& \mathrm{...}& \frac{\&part;N_n}{\&part;x}\\ \frac{\&part;N_1}{\&part;y}& \frac{\&part;N_2}{\&part;y}& \mathrm{...}& \frac{\&part;N_n}{\&part;y}\\ \frac{\&part;N_1}{\&part;z}& \frac{\&part;N_2}{\&part;z}& \mathrm{...}& \frac{\&part;N_n}{\&part;z}\end{array}\right],\hat{T}=\left[\begin{array}{cccc}\frac{\&part;N_1}{\&part;\mathrm{\&xi;}}& \frac{\&part;N_2}{\&part;\mathrm{\&xi;}}& \mathrm{...}& \frac{\&part;N_n}{\&part;\mathrm{\&xi;}}\\ \frac{\&part;N_1}{\&part;\mathrm{\&eta;}}& \frac{\&part;N_2}{\&part;\mathrm{\&eta;}}& \mathrm{...}& \frac{\&part;N_n}{\&part;\mathrm{\&eta;}}\\ \frac{\&part;N_1}{\&part;\mathrm{\&zeta;}}& \frac{\&part;N_2}{\&part;\mathrm{\&zeta;}}& \mathrm{...}& \frac{\&part;N_n}{\&part;\mathrm{\&zeta;}}\end{array}\right]---\left(8\right)$

根据函数求导链式法则得到:

$\hat{T}=J_{\mathrm{\&xi;}}T---\left(9\right)$

根据形函数公式x(ξ)＝P_e ^TN得到

$J_{\mathrm{\&xi;}}=\hat{T}X---\left(10\right)$

P_e是单元控制点坐标矩阵：

$P_e=\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ .& .& .\\ .& .& .\\ .& .& .\\ x_n& y_n& z_n\end{array}\right]---\left(11\right)$

由上面两式求导可得

$\frac{\&part;T}{\&part;{\mathrm{\&alpha;}}_i}=-T\frac{\&part;P_e}{\&part;{\mathrm{\&alpha;}}_i}T---\left(12\right)$

(3)接着需要计算

$\frac{\&part;\left|J\right|}{\&part;{\mathrm{\&alpha;}}_i}=\frac{\&part;\left(\right|J_{\widehat{\mathrm{\&xi;}}}|\&CenterDot;|J_{\widetilde{\mathrm{\&xi;}}}\left|\right)}{\&part;{\mathrm{\&alpha;}}_i}=\left|J_{\widehat{\mathrm{\&xi;}}}\right|\&CenterDot;\frac{\&part;\left(\right|J_{\widehat{\mathrm{\&xi;}}}\left|\right)}{\&part;{\mathrm{\&alpha;}}_i}=\left|J_{\widehat{\mathrm{\&xi;}}}\right|\&CenterDot;tr\left(J^{-1}\hat{T}\frac{\&part;P}{\&part;{\mathrm{\&alpha;}}_i}\right)=\left|J_{\widehat{\mathrm{\&xi;}}}\right|\&CenterDot;tr\left(T\frac{\&part;P_e}{\&part;{\mathrm{\&alpha;}}_i}\right)---\left(13\right)$

其中，tr为矩阵对角线元素之和计算符号，tr(A)代表矩阵A对角线元素之和。

载荷表达为:

$f_e=\underset{{\widetilde{\mathrm{\&Omega;}}}_e}{\&Integral;}{N}^{T}f\left|J\right|d{\widetilde{\mathrm{\&Omega;}}}_e---\left(14\right)$

其中N为形函数矩阵：

$N=\left[\begin{array}{ccccccc}{N}_1& & & \mathrm{...}& N_n& & \\ & N_1& & \mathrm{...}& & N_n& \\ & & N_1& \mathrm{...}& & & N_n\end{array}\right]---\left(15\right)$

则载荷对设计变量导数为：

$\frac{\&part;f_e}{\&part;{\mathrm{\&alpha;}}_i}=\underset{{\widetilde{\mathrm{\&Omega;}}}_e}{\&Integral;}{N}^T\left(\frac{\&part;f}{\&part;{\mathrm{\&alpha;}}_i}\right|J|+f\frac{\&part;\left|J\right|}{\&part;{\mathrm{\&alpha;}}_i})d{\widetilde{\mathrm{\&Omega;}}}_e---\left(16\right)$

在模型所加载荷独立于优化设计参数α_i，故

对平衡方程KU＝F两边同时对设计变量求导得

$\frac{du}{{\mathrm{d\&alpha;}}_i}=-K^{-1}\frac{\&part;K}{\&part;{\mathrm{\&alpha;}}_i}u---\left(17\right)$

其中，B为应变矩阵；D为弹性矩阵；N_i为NURBS基函数表示的形函数；J_ξ为雅克比矩阵；P_e为单元控制点坐标矩阵f_e为载荷。

下面以两个具体实验例来对上述实施例中的一般性的方法进行具体说明，以便更清楚的理解上述方法的设计思路。

实验例一：带脊悬臂梁的形状优化

其形状优化算法包括下述步骤：

步骤1，建立带脊悬臂梁的特征框架模型，如图3(a)和图3(b)所示，该特征框架模型的特征尺寸包括脊的宽度W、脊的高度H、整体厚度t、以及脊的长度L，如图3(c)所示，该模型载荷为Z＝0端面受到均匀沿y轴负向的拉力F＝-50N，边界条件为末端面完全固定；

步骤2，根据式(18)～(21)得到悬臂脊参数化控制点坐标矩阵：

${x_1}^{\&prime;}=x_1+\frac{W-W_0}{W_0}(x_1-x_0)---\left(18\right)$

此处，x₀为脊顶点x坐标，x₁为脊上控制点x坐标，x₁′为更新后的x₁，W₀为顶点脊宽。

脊上控制点y坐标：

满足x＝x₁, (19)

x∈控制点(x,y),y₁∈控制点(x₁,y₁)

此处，y₁为脊底层y坐标，y为脊上控制点y坐标，y′为更新后的y，H₀为底层脊高。

脊外控制点y坐标：

${y_1}^{\&prime;}=y_1+\frac{t-t_0}{t_0}{y}_1---\left(20\right)$

此处，y₁为脊外控制点y坐标，y′为更新后的y₁。

控制点z坐标：

$z^{\&prime;}=z_1+\frac{L-L_0}{L_0}{z}_1---\left(21\right)$

此处，z₁为脊上控制点初始z坐标，z′为更新后的z₁，，L₀为L初始尺寸。

步骤3，根据式(22)得到悬臂脊参数化控制点坐标对特征尺寸的灵敏度矩阵，

$\begin{array}{c}\frac{dP}{dW}=\frac{dP}{{{\mathrm{dx}}_1}^{\&prime;}}\frac{{\mathrm{dx}}_1^{\&prime;}}{dW}=\frac{x_1-x_0}{W_0}\\ \frac{dP}{dH}=\frac{dP}{{\mathrm{dy}}^{\&prime;}}\frac{{\mathrm{dy}}^{\&prime;}}{dH}=\frac{y_1-y_0}{H_0}\\ \frac{dP}{dt}=\frac{dP}{{\mathrm{dy}}^{\&prime;}}\frac{dz}{dt}=\frac{y-y_1}{t_0}\\ \frac{dP}{dL}=\frac{dP}{{\mathrm{dz}}^{\&prime;}}\frac{{\mathrm{dz}}^{\&prime;}}{dL}=\frac{z-z_1}{L_0}\end{array}---\left(22\right)$

步骤4，根据式(23)中所示的模型特征尺寸初始值及约束条件对上述四个特征参数依次进行迭代(迭代历史图如图4所示)，直到满足给定的迭代终止条件。

$\left[\begin{array}{c}{W}_0\\ t_0\\ H_0\\ L_0\end{array}\right]=\left[\begin{array}{c}8\\ 1\\ 3.5\\ 7\end{array}\right],\left[\begin{array}{c}5\\ 0.3\\ 1.2\\ 6\end{array}\right]\&le;\left[\begin{array}{c}W\\ t\\ H\\ L\end{array}\right]\&le;\left[\begin{array}{c}10\\ 1.2\\ 5\\ 10\end{array}\right]---\left(23\right)$

优化前后的结果对比，请见表1：

表1设计目标与约束优化前后对比表

类别	柔度(10-3J)	体积(mm3)	厚度(mm)	脊宽度(mm)	脊高度(mm)	长度(mm)
							优化前	2063	126	1	8	3.5	7
优化后	2114.9	90	0.8	7.8	1.2	6

实验例二：疏浚绞刀臂的形状优化算法

步骤1，建立疏浚绞刀臂的特征框架模型，如图5所示，该特征框架模型的特征尺寸包括绞刀臂下底面轴向宽度bW和轴向长度bL，上顶面轴向宽度tW和轴向长度tL，绞刀臂高度H，载荷与边界条件为：绞刀臂上底面与下底面固定约束，在指定控制点处施加集中力。

步骤2，特征尺寸类型为长度时，映射函数推导。设空间中有任意两点分别为(x₀,y₀,z₀)与(x₁,y₁,z₁)，沿X轴、Y轴、Z轴方向特征尺寸由W₀,L₀,H₀变为W,L,H，将(x₀,y₀,z₀)作为定点，设动点(x₁,y₁,z₁)由于特征尺寸变化得到新的坐标为(x₁',y₁',z₁')，那么映射函数为：

$\begin{array}{c}{y_1}^{\&prime;}=y_1+\frac{W-W_0}{W_0}(y_1-y_0)\\ {x_1}^{\&prime;}=x_1+\frac{L-L_0}{L_0}(x_1-x_0)\\ {z_1}^{\&prime;}=z_1+\frac{H-H_0}{H_0}(z_1-z_0)\end{array}---\left(24\right)$

写成矩阵形式为：

$P=\left[\begin{array}{ccc}{x}_1+\frac{L-L_0}{L_0}(x_1-x_0)& y_1+\frac{W-W_0}{W_0}(y_1-y_0)& z_1+\frac{H-H_0}{H_0}(z_1-z_0)\\ x_2+\frac{L-L_0}{L_0}(x_1-x_0)& y_2+\frac{W-W_0}{W_0}(y_1-y_0)& z_2+\frac{H-H_0}{H_0}(z_1-z_0)\\ .& .& .\\ .& .& .\\ .& .& .\\ x_{n-1}+\frac{L-L_0}{L_0}(x_{n-1}-x_0)& y_{n-1}+\frac{W-W_0}{W_0}(y_{n-1}-y_0)& z_{n-1}+\frac{H-H_0}{H_0}(z_{n-1}-z_0)\\ x_n+\frac{L-L_0}{L_0}(x_n-x_0)& y_n+\frac{W-W_0}{W_0}(y_n-y_0)& z_n+\frac{H-H_0}{H_0}(z_n-z_0)\end{array}\right]---\left(25\right)$

其中，n为单个单元中的控制点总数，根据该式能够得到绞刀臂参数化控制点坐标矩阵。

步骤3，根据式(26)得到疏浚绞刀臂参数化控制点坐标对特征尺寸敏度矩阵，

$\begin{array}{c}\frac{dP}{dW}=\frac{dP}{{{\mathrm{dy}}_1}^{\&prime;}}\frac{dy}{dW}=\frac{y_1-y_0}{W_0}\\ \frac{dP}{dL}=\frac{dP}{{{\mathrm{dx}}_1}^{\&prime;}}\frac{dx}{dL}=\frac{x_1-x_0}{L_0}\\ \frac{dP}{dH}=\frac{dP}{{{\mathrm{dz}}_1}^{\&prime;}}\frac{dz}{dH}=\frac{z_1-z_0}{H_0}\end{array}---\left(26\right)$

步骤4，根据式(27)中所示的模型特征尺寸初始值及约束条件对五个特征参数依次进行迭代(迭代历史图如图6所示)，直到满足给定的迭代终止条件。

$\left[\begin{array}{c}{\mathrm{bW}}_0\\ {\mathrm{bL}}_0\\ {\mathrm{tW}}_0\\ {\mathrm{tL}}_0\\ H_0\end{array}\right]=\left[\begin{array}{c}9.1\\ 24.27\\ 6.8\\ 4.27\\ 117.29\end{array}\right],\left[\begin{array}{c}5.5\\ 19.8\\ 2.4\\ 1.8\\ 112.8\end{array}\right]\&le;\left[\begin{array}{c}bW\\ bL\\ tW\\ tL\\ H\end{array}\right]\&le;\left[\begin{array}{c}13.5\\ 29.6\\ 10.2\\ 7.6\\ 122.79\end{array}\right]---\left(27\right)$

优化前后的结果对比，请见表2：

表2设计目标与约束优化前后对比表

类别	柔度(10-3J)	体积(mm3)	bL(mm)	bW(mm)	tL(mm)	tW(mm)	H(mm)
								原始模型	578	18628	24.27	9.1	4.27	6.8	117.29
优化模型	520	18200	29.6	6.3	7.6	6.8	112.8

通过对比两个例子的优化前后变化，在模型体积减小的情况下，模型刚度也得到了一定程度的提高，而且优化效率较高，使得产品性能得到改善。

Claims (9)

1.一种三维产品模型的形状优化算法，用于对三维产品模型的形状进行优化设计，其特征在于，包括以下步骤：

(一)建立特征框架模型

采用反求方法或者尺寸测量方法，构建包含产品模型点、线、面特征的自由形体特征框架模型；选取部分特征参数作为优化设计参数，并根据这些优化设计参数和其他一些已知参数，构建出所述自由形体特征框架模型的特征曲线；给所述特征曲线添加约束条件得到参数化的特征框架，通过特征框架尺寸驱动产品模型的变形，

(二)建立分析模型

采用一定方法获取所述特征框架模型的表面控制点和实体内部控制点，并建立优化设计参数和控制点之间的关联，

(三)建立控制点坐标对优化设计参数的灵敏度矩阵

给所述特征框架模型施加边界条件或约束条件，利用等几何分析的方法得到每个控制点的物性参数；利用所述控制点和优化设计参数的关系，从而构建出优化设计参数和物性参数的关系；根据该关系求取物性参数针对于优化设计参数的导数关系，即得到灵敏度矩阵；对所述灵敏度矩阵进行求解，确定优化算法，

(四)利用选用的优化算法，依次对每个优化设计参数和对应的特征框架模型进行更新，直到满足给定的迭代终止条件。

2.根据权利要求1所述三维产品模型的形状优化算法，其特征在于：

其中，所述一定方法为离散孔斯插值法、调和函数法、凸组合插值法中的任意一种。

3.根据权利要求1所述三维产品模型的形状优化算法，其特征在于：

其中，每个控制点的物性参数为所述控制点的位移和受力。

4.根据权利要求1所述三维产品模型的形状优化算法，其特征在于：

其中，所述三维产品模型的优化目标为模型的刚度最大，所述特征框架模型施的约束条件为模型体积小于特定值，

所述三维产品模型的结构柔度对优化设计参数的导数为：

$\frac{d\left(f^{T}u\right)}{{\mathrm{d\&alpha;}}_i}=f^T\frac{\&part;u}{\&part;{\mathrm{\&alpha;}}_i}=-f^T{K}^{-1}\frac{\&part;K}{\&part;{\mathrm{\&alpha;}}_i}u=-u^T\frac{\&part;K}{\&part;{\mathrm{\&alpha;}}_i}u---\left(1\right)$

所述三维产品模型的体积或面积对优化设计参数的导数为：

$\frac{ds}{{\mathrm{d\&alpha;}}_i}=\underset{e=1}{\overset{n_{el}}{\&Sigma;}}\frac{d\left|J\right|}{{\mathrm{d\&alpha;}}_i}d\widetilde{\mathrm{\&Omega;}}---\left(2\right)$

式中，f为载荷列阵；u为位移列阵；α_i为优化设计参数；K为三维产品模型的整体刚度；S为三维产品模型的体积或面积；|J|为雅克比行列式；Ω为分析域。

5.根据权利要求4所述三维产品模型的形状优化算法，其特征在于：

其中，所述灵敏度矩阵的求解方法包括：

(1)整体刚度矩阵K是由单元刚度矩阵K_e组装而来，因此转化为求单元刚度矩阵K_e对优化设计参数α_i的导数问题：

$\frac{\&part;K_e}{\&part;{\mathrm{\&alpha;}}_i}=\underset{{\widetilde{\mathrm{\&Omega;}}}_e}{\&Integral;}\frac{\&part;B^T}{\&part;\mathrm{\&alpha;}}DB\left|J\right|d{\widetilde{\mathrm{\&Omega;}}}_e+\underset{{\widetilde{\mathrm{\&Omega;}}}_e}{\&Integral;}{B}^{T}D\frac{\&part;B}{\&part;\mathrm{\&alpha;}}\left|J\right|d{\widetilde{\mathrm{\&Omega;}}}_e+\underset{{\widetilde{\mathrm{\&Omega;}}}_e}{\&Integral;}{B}^{T}DB\frac{\&part;\left|J\right|}{\&part;\mathrm{\&alpha;}}d{\widetilde{\mathrm{\&Omega;}}}_e---\left(3\right)$

(2)首先计算B对α_i的导数：

B＝[B₁ B₂ ... B_i ... B_n] (4)

$B_i=\left[\begin{array}{ccc}\frac{\&part;N_i}{\&part;x}& 0& 0\\ 0& \frac{\&part;N_i}{\&part;y}& 0\\ 0& 0& \frac{\&part;N_i}{\&part;z}\\ \frac{\&part;N_i}{\&part;y}& \frac{\&part;N_i}{\&part;x}& 0\\ 0& \frac{\&part;N_i}{\&part;z}& \frac{\&part;N_i}{\&part;y}\\ \frac{\&part;N_i}{\&part;z}& 0& \frac{\&part;N_i}{\&part;x}\end{array}\right]---\left(5\right)$

定义两个矩阵

$T=\left[\begin{array}{cccc}\frac{\&part;N_1}{\&part;x}& \frac{\&part;N_2}{\&part;x}& \mathrm{...}& \frac{\&part;N_n}{\&part;x}\\ \frac{\&part;N_1}{\&part;y}& \frac{\&part;N_2}{\&part;y}& \mathrm{...}& \frac{\&part;N_n}{\&part;y}\\ \frac{\&part;N_1}{\&part;z}& \frac{\&part;N_2}{\&part;z}& \mathrm{...}& \frac{\&part;N_n}{\&part;z}\end{array}\right],\hat{T}=\left[\begin{array}{cccc}\frac{\&part;N_1}{\&part;\mathrm{\&xi;}}& \frac{\&part;N_1}{\&part;\mathrm{\&xi;}}& \mathrm{...}& \frac{\&part;N_n}{\&part;\mathrm{\&xi;}}\\ \frac{\&part;N_1}{\&part;\mathrm{\&eta;}}& \frac{\&part;N_2}{\&part;\mathrm{\&eta;}}& \mathrm{...}& \frac{\&part;N_n}{\&part;\mathrm{\&eta;}}\\ \frac{\&part;N_1}{\&part;\mathrm{\&zeta;}}& \frac{\&part;N_2}{\&part;\mathrm{\&zeta;}}& \mathrm{...}& \frac{\&part;N_n}{\&part;\mathrm{\&zeta;}}\end{array}\right]---\left(6\right)$

根据函数求导链式法则得到:

$\hat{T}=J_{\mathrm{\&xi;}}T---\left(7\right)$

根据形函数公式得到

$J_{\mathrm{\&xi;}}=\hat{T}{P}_e---\left(8\right)$

P_e是单元控制点坐标矩阵：

$P_e=\left[\begin{array}{ccc}{x}_1& y_1& z_1\\ .& .& .\\ .& .& .\\ .& .& .\\ x_n& y_n& z_n\end{array}\right]---\left(9\right)$

由上面两式求导可得

$\frac{\&part;T}{\&part;{\mathrm{\&alpha;}}_i}=-T\frac{\&part;P_e}{\&part;{\mathrm{\&alpha;}}_i}T---\left(10\right)$

(3)接着需要计算

$\frac{\&part;\left|J\right|}{\&part;{\mathrm{\&alpha;}}_i}=\frac{\&part;\left(\right|J_{\widehat{\mathrm{\&xi;}}}|\&CenterDot;|J_{\widetilde{\mathrm{\&xi;}}}\left|\right)}{\&part;{\mathrm{\&alpha;}}_i}=\left|J_{\widehat{\mathrm{\&xi;}}}\right|\&CenterDot;\frac{\&part;\left(\right|J_{\widehat{\mathrm{\&xi;}}}\left|\right)}{\&part;{\mathrm{\&alpha;}}_i}=\left|J_{\widehat{\mathrm{\&xi;}}}\right|\&CenterDot;tr\left(J^{-1}\hat{T}\frac{\&part;P_e}{\&part;{\mathrm{\&alpha;}}_i}\right)=\left|J_{\widehat{\mathrm{\&xi;}}}\right|\&CenterDot;tr\left(T\frac{\&part;P_e}{\&part;{\mathrm{\&alpha;}}_i}\right)---\left(11\right)$

其中，tr为矩阵对角线元素之和计算符号，tr(A)代表矩阵A对角线元素之和。

载荷表达为:

$f_e=\underset{{\widetilde{\mathrm{\&Omega;}}}_e}{\&Integral;}{N}^{T}f\left|J\right|d{\widetilde{\mathrm{\&Omega;}}}_e---\left(12\right)$

其中N为形函数矩阵：

$N=\left[\begin{array}{ccccccc}{N}_1& & & ...& N_n& & \\ & N_1& & ...& & N_n& \\ & & N_1& ...& & & N_n\end{array}\right]---\left(13\right)$

则载荷对设计变量导数为：

$\frac{\&part;f_e}{\&part;{\mathrm{\&alpha;}}_i}=\underset{{\widetilde{\mathrm{\&Omega;}}}_e}{\&Integral;}{N}^T\left(\frac{\&part;f}{\&part;{\mathrm{\&alpha;}}_i}\right|J|+f\frac{\&part;\left|J\right|}{\&part;{\mathrm{\&alpha;}}_i})d{\widetilde{\mathrm{\&Omega;}}}_e---\left(14\right)$

在模型所加载荷独立于优化设计参数α_i，故

对平衡方程KU＝F两边同时对设计变量求导得

$\frac{du}{{\mathrm{d\&alpha;}}_i}=-K^{-1}\frac{\&part;K}{\&part;{\mathrm{\&alpha;}}_i}u---\left(15\right)$

其中，B为应变矩阵；D为弹性矩阵；N_i为NURBS基函数表示的形函数；J_ξ为雅克比矩阵；P_e为单元控制点坐标矩阵f_e为载荷。

6.权利要求1～5中任一项所述的三维产品模型的形状优化算法在带脊悬臂梁形状优化中的用途。

7.根据权利要求6所述的带脊悬臂梁的形状优化算法，其特征在于，包括以下步骤：

步骤1，建立带脊悬臂梁的特征框架模型，该特征框架模型的特征尺寸包括脊的宽度W、脊的高度H、整体厚度t、以及脊的长度L，该模型载荷为Z＝0端面受到均匀沿y轴负向的拉力F，边界条件为末端面完全固定；

步骤2，根据式(16)～(19)得到悬臂脊参数化控制点坐标矩阵：

${x_1}^{\&prime;}=x_1+\frac{W-W_0}{W_0}(x_1-x_0)---\left(16\right)$

此处，x₀为脊顶点x坐标，x₁为脊上控制点x坐标，x₁′为更新后的x₁，W₀为顶点脊宽，

脊上控制点y坐标：

满足x＝x₁,(17)

x∈控制点(x,y),y₁∈控制点(x₁,y₁)

此处，y₁为脊底层y坐标，y为脊上控制点y坐标，y′为更新后的y，

H₀为底层脊高，

脊外控制点y坐标：

${y_1}^{\&prime;}=y_1+\frac{t-t_0}{t_0}{y}_1---\left(18\right)$

此处，y₁为脊外控制点y坐标，y′为更新后的y₁，

控制点z坐标：

$z^{\&prime;}=z_1+\frac{L-L_0}{L_0}{z}_1---\left(19\right)$

此处，z₁为脊上控制点初始z坐标，z′为更新后的z₁，，L₀为L初始尺寸。

步骤3，根据式(20)得到悬臂脊参数化控制点坐标对特征尺寸的灵敏度矩阵，

$\begin{array}{c}\frac{dP}{dW}=\frac{dP}{{{\mathrm{dx}}_1}^{\&prime;}}\frac{{\mathrm{dx}}_1^{\&prime;}}{dW}=\frac{x_1-x_0}{W_0}\\ \frac{dP}{dH}=\frac{dP}{{\mathrm{dy}}^{\&prime;}}\frac{{\mathrm{dy}}^{\&prime;}}{dH}=\frac{y_1-y_0}{H_0}\\ \frac{dP}{dt}=\frac{dP}{{\mathrm{dy}}^{\&prime;}}\frac{dz}{dt}=\frac{y-y_1}{t_0}\\ \frac{dP}{dL}=\frac{dP}{{\mathrm{dz}}^{\&prime;}}\frac{{\mathrm{dz}}^{\&prime;}}{dL}=\frac{z-z_0}{L_0}\end{array}---\left(20\right)$

步骤4，根据式(21)中所示的模型特征尺寸初始值及约束条件对柔度以及上述四个特征参数依次进行迭代，直到满足给定的迭代终止条件。

$\left[\begin{array}{c}{W}_0\\ t_0\\ H_0\\ L_0\end{array}\right]=\left[\begin{array}{c}8\\ 1\\ 3.5\\ 7\end{array}\right],\left[\begin{array}{c}5\\ 0.3\\ 1.2\\ 6\end{array}\right]\&le;\left[\begin{array}{c}W\\ t\\ H\\ L\end{array}\right]\&le;\left[\begin{array}{c}10\\ 1.2\\ 5\\ 10\end{array}\right]---\left(21\right)$

8.权利要求1～5中任一项所述的三维产品模型的形状优化算法在疏浚绞刀臂形状优化中的用途。

9.根据权利要求8所述的疏浚绞刀臂的形状优化算法，其特征在于，包括以下步骤：

步骤1，建立疏浚绞刀臂的特征框架模型，该特征框架模型的特征尺寸包括绞刀臂下底面轴向宽度bW和轴向长度bL，上顶面轴向宽度tW和轴向长度tL，绞刀臂高度H，载荷与边界条件为：绞刀臂上底面与下底面固定约束，在指定控制点处施加集中力；

步骤2，根据式(22)得到疏浚绞刀臂参数化控制点坐标矩阵：

$P=\left[\begin{array}{ccc}{x}_1+\frac{L-L_0}{L_0}(x_1-x_0)& y_1+\frac{W-W_0}{W_0}(y_1-y_0)& z_1+\frac{H-H_0}{H_0}(z_1-z_0)\\ x_2+\frac{L-L_0}{L_0}(x_1-x_0)& y_2+\frac{W-W_0}{W_0}(y_1-y_0)& z_2+\frac{H-H_0}{H_0}(z_1-z_0)\\ .& .& .\\ .& .& .\\ .& .& .\\ x_{n-1}+\frac{L-L_0}{L_0}(x_{n-1}-x_0)& y_{n-1}+\frac{W-W_0}{W_0}(y_{n-1}-y_0)& z_{n-1}+\frac{H-H_0}{H_0}(z_{n-1}-z_0)\\ x_n+\frac{L-L_0}{L_0}(x_n-x_0)& y_n+\frac{W-W_0}{W_0}(y_n-y_0)& z_n+\frac{H-H_0}{H_0}(z_n-z_0)\end{array}\right]---\left(22\right)$

此处，(x₀,y₀,z₀)为空间中任意设定的定点，(x₁,y₁,z₁)为动点，W₀,L₀,H₀为与(x₀,y₀,z₀)相对应地沿X轴、Y轴、Z轴方向特征尺寸，W,L,H为与(x₁,y₁,z₁)相对应地沿X轴、Y轴、Z轴方向特征尺寸，

步骤3，根据式(23)得到疏浚绞刀臂参数化控制点坐标对特征尺寸敏度矩阵，

$\begin{array}{c}\frac{dP}{dW}=\frac{dP}{{{\mathrm{dy}}_1}^{\&prime;}}\frac{dy}{dW}=\frac{y_1-y_0}{W_0}\\ \frac{dP}{dL}=\frac{dP}{{{\mathrm{dx}}_1}^{\&prime;}}\frac{dx}{dL}=\frac{x_1-x_0}{L_0}\\ \frac{dP}{dH}=\frac{dP}{{{\mathrm{dz}}_1}^{\&prime;}}\frac{dz}{dH}=\frac{z_1-z_0}{H_0}\end{array}---\left(23\right)$

步骤4，根据式(24)中所示的模型特征尺寸初始值及约束条件对五个特征参数依次进行迭代，直到满足给定的迭代终止条件。

$\left[\begin{array}{c}{\mathrm{bW}}_0\\ {\mathrm{bL}}_0\\ {\mathrm{tW}}_0\\ {\mathrm{tL}}_0\\ H_0\end{array}\right]=\left[\begin{array}{c}9.1\\ 24.27\\ 6.8\\ 4.27\\ 117.29\end{array}\right],\left[\begin{array}{c}5.5\\ 19.8\\ 2.4\\ 1.8\\ 112.8\end{array}\right]\&le;\left[\begin{array}{c}bW\\ bL\\ tW\\ tL\\ H\end{array}\right]\&le;\left[\begin{array}{c}13.5\\ 29.6\\ 10.2\\ 7.6\\ 122.79\end{array}\right]---\left(24\right)$