CN103970952A - 一种面向三维制造技术的快速应力分析方法 - Google Patents

Abstract

本发明涉及一种面向三维制造技术的快速应力分析方法，主要包括以下步骤：对输入网格模型的预处理；Boundary mode近似参数函数曲面系数的预计算；参数化stiffness矩阵的预计算；在用户进行编辑时即时更新应力求解系统的运算数据并实时估算应力分布情况；编辑完成后提供应力分布情况的准确计算。发明有益的效果是：本发明通过将网格模型的编辑和快速应力估算技术相结合，采用参数化、逼近和解空间降维等方法，发明了一种可以实时预览编辑操作对模型应力分布所造成的影响的技术，为模型设计过程提供直观的辅助，缩短了三维模型打印技术的设计和制造周期。

Description

一种面向三维制造技术的快速应力分析方法

技术领域

本发明涉及3D打印制造领域，尤其是一种面向三维制造技术的快速应力分析方法。

背景技术

3D打印是当下热门的研究领域，其研究工作主要集中在打印结果的物理形变(BICKEL,B.,BACHER,M.,OTADUY,M.A.,LEE,H.R.,PFISTER,H.,GROSS,M.,AND MATUSIK,W.2010.Design and fabrication of materials with desired deformation behavior.ACM Trans.Graph.29,4(July),63:1–63:10；SKOURAS,M.,THOMASZEWSKI,B.,COROS,S.,BICKEL,B.,AND GROSS,M.2013.Computational design of actuated deformable characters.ACMTrans.Graph.32,4(July),82:1–82:10.),打印结果的关节结构(CAL`I，J.,CALIAN,D.A.,AMATI,C.,KLEINBERGER,R.,STEED,A.,KAUTZ,J.,AND WEYRICH,T.2012.3d-printing of non-assembly,articulated models.ACM Trans.Graph.31,6(Nov.),130:1–130:8；B¨A CHER,M.,BICKEL,B.,JAMES,D.L.,AND PFISTER,H.2012.Fabricating articulated characters from skinned meshes.ACM Trans.Graph.31,4(July),47:1–47:9)，打印结果的机械活动性(ZHU,L.,XU,W.,SNYDER,J.,LIU,Y.,WANG,G.,AND GUO,B.2012.Motion-guided mechanical toy modeling.ACM Trans.Graph.31,6(Nov.),127:1–127:10；COROS,S.,THOMASZEWSKI,B.,NORIS,G.,SUEDA,S.,FORBERG,M.,SUMNER,R.W.,MATUSIK,W.,AND BICKEL,B.2013.Computational design ofmechanical characters.ACM Trans.Graph.32,4(July),83:1–83:12；CEYLAN,D.,LI,W.,MITRA,N.J.,AGRAWALA,M.,AND PAULY,M.2013.Designing and fabricatingmechanical automata from mocap sequences.ACM Trans.Graph.32,6(Nov.),186:1–186:11.)，以及打印结果的外观(DONG,Y.,WANG,J.,PELLACINI,F.,TONG,X.,AND GUO,B.2010.Fabricating spatially-varying subsurface scattering.ACM Trans.Graph.29,4(July),62:1–62:10；LAN,Y.,DONG,Y.,PELLACINI,F.,AND TONG,X.2013.Bi-scale appearance fabrication.ACM Trans.Graph.32,4(July),145:1–145:12；CHEN,D.,LEVIN,D.I.W.,DIDYK,P.,SITTHI-AMORN,P.,AND MATUSIK,W.2013.Spec2fab:A reducer-tuner model for translating specifications to3d prints.ACM Trans.Graph.32,4(July),135:1–135:10)。

制造相关设计的研究重点是简化制造过程,(LIU,Y.,POTTMANN,H.,WALLNER,J.,YANG,Y.-L.,AND WANG,W.2006.Geometric modeling with conical meshes and developablesurfaces.ACM Trans.Graph.25,3(July),681–689；LIU,Y.,XU,W.,WANG,J.,ZHU,L.,GUO,B.,CHEN,F.,AND WANG,G.2011.General planar quadrilateral mesh designusing conjugate direction field.ACM Trans.Graph.30,6(Dec.),140:1–140:10)开发了平面四边形网格算法以减少建筑结构表面建造的开销。(MCCRAE,J.,SINGH,K.,ANDMITRA,N.J.2011.Slices:A shape-proxy based on planar sections.ACM Trans.Graph.30,6(Dec.),168:1–168:12.；HILDEBRAND,K.,BICKEL,B.,AND ALEXA,M.2012.Crdbrd:Shape fabrication by sliding planar slices.Comp.Graph.Forum31,2pt3(May),583–592)开发了将3D形状转变为平面切片图形的算法。(LAU,M.,OHGAWARA,A.,MITANI,J.,AND IGARASHI,T.2011.Converting3d furniture models to fabricatable parts andconnectors.ACM Trans.Graph.30,4(July),85:1–85:6)开发了将3D家具模型按照制造手段的约束而分拆成部分的方法。(MORI,Y.,AND IGARASHI,T.2007.Plushie:Aninteractive design system for plush toys.ACM Trans.Graph.26,3(July))和(UMETANI,N.,KAUFMAN,D.M.,IGARASHI,T.,AND GRINSPUN,E.2011.Sensitive couture forinteractive garment modeling and editing.ACM Trans.Graph.30,4(July),90:1–90:12)分别在玩具制造和服装设计上将快速物理模拟技术集成到设计系统中，使用户可以预测所设计物品的物理属性。(UMETANI,N.,IGARASHI,T.,AND MITRA,N.J.2012.Guided exploration of physically valid shapes for furniture design.ACM Trans.Graph.31,4(July),86:1–86:11)在家具设计系统中集成了物体受力薄弱点分析的技术。

区域划分有限元算法：(TOSELLI,A.,AND WIDLUND,O.2004.Domain DecompositionMethods.Springer)根据分而治之的思想，将区域划分用在偏微分方法的求解上。(MANDEL,J.,AND TEZAUR,R.2001.On the convergence of a dual-primal substructuring method.Numerische Mathematik88,3,543–558；FARHAT,C.,AND ROUX,F.-X.1991.A methodof finite element tearing and interconnecting and its parallel solution algorithm.International Journal for Numerical Methods in Engineering32,6,1205–1227；FARHAT,C.,LESOINNE,M.,LETALLEC,666P.,PIERSON,K.,AND RIXEN,D.2001.Feti-dp:adual-primal unified feti method part i:A faster alternative to the two-level fetimethod.International Journal for Numerical Methods in Engineering50,7,1523–1544.)将此技术应用在应力分布分析中。(BARBIC,J.,AND ZHAO,Y.2011.Real-time large-deformation substructuring.In ACM SIGGRAPH2011Papers,SIGGRAPH’11,91:1–91:8,KIM,T.,AND JAMES,D.L.2011.Physics-based character skinning usingmulti-domain subspace deformations.In Proceedings of the2011ACMSIGGRAPH/Eurographics Symposium on Computer Animation,SCA’11,63–72，YANG,Y.,XU,W.,GUO,X.,ZHOU,K.,AND GUO,B.2013.Boundary-aware multidomain subspacedeformation.Visualization and Computer Graphics,IEEE Transactions on19,10,1633–1645)致力于简化区域划分技术以换取求解速度上的提升。

发明内容

本发明提供一种面向三维制造技术的快速应力分析方法，这是一种将网格模型编辑和应力分析相结合的方法，即能对应力分布情况进行实时预测，又能对应力分布情况进行准确求解，从而辅助用户设计安全可靠的3D打印模型。

本发明解决其技术问题采用的技术方案：一种面向三维制造技术的快速应力分析方法，主要包括以下步骤：

1)对输入网格模型的预处理；

2)Boundary mode近似参数函数曲面系数的预计算；

3)参数化stiffness矩阵的预计算；

4)在用户进行编辑时即时更新应力求解系统的运算数据并实时估算应力分布情况；

5)编辑完成后提供应力分布情况的准确计算。

作为优选，步骤2)中提出了一种Boundary mode的快速估算方法，通过将Boundary mode的元素视作高维曲面，用参数曲面函数通过最小二乘的方式逼近该曲面。将逼近得到的参数曲面系数保存成矩阵的格式，当网格模型的缩放系数发生变化时，可以通过矩阵相加的方式快速得到Boundary mode的估算值。

作为优选，步骤3)中提出了一种将stiffness矩阵参数化的方法，通过将矩阵元素表示成多项式的形式，并对多项式的系数进行预计算。当网格的缩放参数和物理参数发生变化时，便可以快速地得到网格的stiffness矩阵，而不用重新遍历网格模型进行计算。

作为优选，步骤4)中提出了一种基于区域划分的网格模型编辑方法，通过将网格模型划分为子区域，提供针对网格区域进行的编辑功能，从而实现对模型整体进行编辑的功能。

作为优选，步骤4)中提出了一种将网格模型编辑和快速应力近似逼近求解技术相结合的模型设计方法，使得系统用户可以在对网格模型进行编辑的同时，可实时地预览编辑操作对模型造成的应力分布情况的改变，从而辅助用户设计可安全打印的三维网格模型。

作为优选，步骤4)中提出了一种将FETI-DP方法的解空间降维到Boundary mode子空间的快速求解方法。通过将FETI-DP中每个网格域的自由度个数降低到边界顶点的个数，从而降低了整个系统自由度的数目，使得系统可以快速地得到原系统的近似解。

发明有益的效果是：本发明通过将网格模型的编辑和快速应力估算技术相结合，采用参数化、逼近和解空间降维等方法，发明了一种可以实时预览编辑操作对模型应力分布所造成的影响的技术，为模型设计过程提供直观的辅助，缩短了三维模型打印技术的设计和制造周期。

附图说明

图1是参数曲面函数逼近boundary mode高维曲面示意图；

图2是骨骼和区域划分示意图；

图3是编辑前后形变以及应力分布对比示意图；

具体实施方式

下面结合附图对本发明作进一步说明：

一种面向三维制造技术的快速应力分析方法，分为如下几个主要步骤：对输入网格模型的预处理；对Boundary mode近似参数曲面的预计算；参数化stiffness矩阵的预计算；在用户进行编辑时即时更新应力求解系统运算数据并提供应力分布的估算值；编辑完成后提供准确的应力分布情况。

1)网格模型的预处理：本发明所需的输入是一个四面体网格模型，系统从输入的网格模型中自动或人工地抽取出模型的骨骼结构，然后基于所抽取得到的骨骼结构，将网格模型划分为区域，并计算可形变区域顶点的混合权重。

1.1)从四面体模型中抽取骨骼结构：可以采用(BARAN,I.,AND POPOVI′C,J.2007.Automatic rigging and animation of3d characters.ACM Trans.Graph.26,3(July))中的方法，自动从四面体网格模型中抽取骨骼结构，然后视场景需要手工对抽取得到的骨骼结构进行修改，也可以直接在建模软件中手工生成骨骼结构。

1.2)根据骨骼结构将四面体网格划分为区域：如果采用1.1)中的自动方法，结果中已经包含了网格顶点的动画权重。将动画权重中相对于某一骨骼权重为1.0的顶点分别划归同一区域，便可得到绑定于该骨骼的刚体运动网格区域。网格中剩余顶点作为可形变区域。用户也可以手工地自行划分区域，或者对现有区域划分进行修改(如图2所示)。

1.3)根据骨骼和区域赋予可形变区域的顶点混合权重：可形变区域顶点的混合权重可由求解如下系统得到：

Δwⁱ+H(pⁱ-wⁱ)＝0

其中Δ为四面体网格Laplacian算子，wⁱ表示各顶点处相对于骨骼i的混合权重，为该系统所要求解的值。pⁱ为顶点权重的目标值，当顶点位于区域间的邻接面，且该邻接面位于骨骼i上时，将p值设为，请他情况将p设为0。H表示一个对角矩阵，其行数列数为区域内顶点的自由度数目，其对角元素为p值和w值的差。

2)Boundary mode近似参数函数曲面系数的预计算：本发明将Boundary Mode保存成矩阵的形式。Boundary mode中的每个元素可视作一个高维曲面，其坐标轴为三个缩放轴。本发明用参数曲面函数近似逼近原高维曲面(如图1所示)，该函数的公式为：

$\begin{array}{c}{\mathrm{\Φ}}^{\⟨s_x,s_y,s_z\⟩}=\mathrm{\Φ}+c_0{(s_x-1)}^2+c_1{(s_y-1)}^2+c_2{(s_z-1)}^2+c_3(s_x-1)(s_y-1)+\\ c_4(s_x-1)(s_z-1)+c_5(s_y-1)(s_z-1)+c_6(s_x-1)+c_7(s_y-1)+c_8(s_z-1)\end{array}$

该参数曲面函数的参数为s_x,s_y,s_z，即三个缩放系数。Φ为该元素在没有经过编辑(即s_x＝1,s_y＝1,s_z＝1)处的值，c_k为所需进行预计算的参数曲面函数的系数值。

本发明在s_x，s_y,s_z三个坐标轴上各选取一些值作为采样点，采用最小二乘的方法计算出系数c_k。本发明以矩阵的形式保存未经编辑时的Boundary mode值(Φ)，以及计算得到的函数系数(c_k)。提供编辑后的三个缩放系数，计算(s_x-1)等项，通过矩阵累加的方式，即可得到Boundary mode的估算值。

Boundary mode估算装置包括：曲面函数系数保存装置，以矩阵的形式保存预计算得到的参数曲面函数的系数值；Boundary mode估算装置，给定一组缩放系数，以矩阵累加的方式得到Boundary mode的估算值。

3)参数化stiffness矩阵的预计算：本发明将stiffness矩阵中的每一项表示为多项式的形式，通过预计算多项式系数来加速stiffness矩阵的构造。

3.1)矩阵项多项式公式：对于单个四面体元，其stiffness矩阵的公式为：

$K_e=\∫\∫\∫B^T\mathrm{DBdxdydz}$

其中D为根据材料物理系数得到的矩阵，B为将顶点位移转换为strain的矩阵。将该式展开，则矩阵的每一项可以表示为：

$c(i,j,{\mathrm{\θ}}_1,{\mathrm{\θ}}_2)=\∫\frac{{\∂N}_i}{{\∂\mathrm{\θ}}_1}\frac{{\∂N}_j}{{\∂\mathrm{\θ}}_2}\mathrm{dxdydz}$

其中i和j为元素在矩阵中的位置，θ₁,θ₂根据元素位置可分别取符号x,y,z，N为四面体元的形状函数。顶点的坐标可以表示为

x_i＝x'_p+s_x(x′_i-x'_p)

y_i＝y'_p+s_y(y′_i-y'_p)

z_i＝z'_p+s_z(z′_i-z'_p)

式中x′_i表示顶点i的初始坐标值，x'_p表示缩放中心点的坐标值，x_i为缩放后的顶点坐标值。讲上式带入stiffness矩阵项的公式的形状函数中，并进一步展开，最终可将矩阵项表示为：

$\begin{array}{c}c(i,j)=c_0{s}_y^2{s}_z^2+c_1{s}_z^2{s}_x^2+c_2{s}_x^2{s}_y^2+c_3{s}_x^2{s}_y{s}_z+c_4{s}_y^2{s}_z{s}_x+\\ c_5{\mathrm{vs}}_z^2{s}_x{s}_y+c_6{\mathrm{vs}}_y^2{s}_z^2+c_7{\mathrm{vs}}_z^2{s}_x^2+c_8{\mathrm{vs}}_x^2{s}_y^2+c_9{\mathrm{vs}}_x^2{s}_y{s}_z+c_{10}{\mathrm{vs}}_y^2{s}_z{s}_x+c_{11}{\mathrm{vs}}_z^2{s}_x{s}_y\end{array}$

其中s_x,s_y,s_z,为三个缩放系数，v为材料的泊松系数。当四面体网格的初始顶点坐标确定时，通过积分计算可以得到系数c_k。

在三维情况下，采用二次元进行计算的单个四面体元的stiffness矩阵为一个30*30的矩阵。整个四面体网格的stiffness矩阵是一个N*N的大稀疏矩阵，N为四面体网格总自由度的个数。根据单个四面体元顶点的索引号，将每个四面体元的stiffness矩阵分散并累加到包含所有四面体元，顶点统一排序的全局大矩阵中，便完成了整个四面体网格的stiffness矩阵的构建。在本发明将矩阵项表示为多项式的情况下，在累加过程中，将相同参数的系数进行相加，便能够得到整个四面体网格的stiffness矩阵的多项式系数。

3.2)用stiffness矩阵的参数化表示快速生成矩阵：当四面体网格在编辑中被缩放时，本发明只需将缩放后的缩放系数s_x,s_y,s_z,带入到上面推导得出的公式中，通过简单的乘法和加法运算，便可以得到最终的stiffness矩阵，无需再遍历每一个四面体元并进行积分计算，从而提高了stiffness矩阵的构建速度。

通过将stiffness矩阵参数化，快速计算stiffness矩阵的装置包括：多项式系数保存装置，以数组的形式保存1中每个矩阵项的多项式系数；stiffness矩阵构建装置，给定一组缩放参数及物理系数，该装置计算stiffness矩阵每一项的多项式，并将结果保存在一个矩阵中。

4)在用户进行编辑时即时更新应力求解系统的运算数据并实时估算应力分布情况：

4.1)编辑操作：本发明的编辑操作针对区域划分得到的网格域进行。对于刚体运动网格域，能够进行平移、旋转和缩放；对于可形变网格域，能够使用所邻接的刚体网格域中的顶点坐标和混合权重，对可形变网格域的顶点坐标进行混合重建。在这些基本编辑操作的基础上，本发明通过对骨骼结构使用IK技术，提供用户更方便快捷友好的操作界面。

以平移，旋转，缩放，以及顶点坐标混合重建作为基本编辑方法的用户编辑装置包括：以IK方法为基础的骨骼动画编辑装置，用户拖动骨骼末端的控制点，移动从骨骼跟节点至骨骼末端所途径的所有网格域，装置将对其中的刚体形变域施加平移及旋转编辑操作，对可形变域进行顶点坐标混合重建编辑操作；骨骼旋转装置，用户旋转某一骨骼，装置将对该骨骼所影响的所有刚体域施加旋转编辑操作，并对受影响的可形变域施加顶点坐标混合重建操作；刚体运动域缩放装置，用户通过该装置对刚体运动域进行缩放，装置对被缩放的刚体运动域施加缩放编辑操作，对相邻接的受影响可形变域施加顶点坐标混合重建操作。

4.2)即时更新应力求解系统的运算数据：当用户对某一刚体运动网格域进行平移操作时，应力求解系统不需要进行更新；当用户对某一刚体运动网格域进行旋转操作时，本发明更新应力求解系统中该网格域所对应的乘子矩阵；当用户对某一网格域进行缩放操作时，本发明更新应力求解系统中该网格域所对应的stiffness矩阵(通过参数化的方法)，和Boundary mode(通过Boundary mode的近似估算)。

当一个可形变网格域进行顶点坐标混合重建时，本发明重新计算该网格域stiffness矩阵参数化形式的系数，然后重新计算该网格域的stiffness矩阵。可形变网格域通常顶点个数较少，因此对其参数化表示中系数的重新计算不会对系统的整体运算效率产生较大的影响。

在用户编辑过程中即时更新应力求解系统中运算数据的装置包括：域运算数据保存装置，以矩阵和向量的形式保存每个域的求解运算数据；域编辑标记装置，当用户对网格模型进行编辑时，标记编辑所影响的刚体运动域，以及因刚体运动域被编辑而受到影响的可形变域；域系统数据更新装置，按照域编辑标记装置中所标记的数据，更新域运算数据保存装置中所保存的运算数据。

4.3)应力分布情况的快速估算：

4.3.1)应力分布情况的快速求解：本发明对FETIDP方法(FARHAT,C.,LESOINNE,M.,LETALLEC,666P.,PIERSON,K.,AND RIXEN,D.2001.Feti-dp:adual-primal unified feti method part i:A faster alternative to thetwo-level feti method.International Journal for Numerical Methodsin Engineering50,7,1523–1544)进行改进，开发了一种能够快速估算应力分布情况的求解系统。

本方法的改进思想是，用Boundary mode对原系统进行降维。某个网格域d所对应的Boundary mode的形式为：

${\mathrm{\Φ}}^d=\left(\begin{array}{cc}{\mathrm{\Φ}}_{\mathrm{ij}}^d& {\mathrm{\Φ}}_{\mathrm{ic}}^d\\ I_{\mathrm{ff}}& 0\\ 0& I_{\mathrm{cc}}\end{array}\right)$

将stiffness矩阵按照内部顶点，边界顶点和corner point的分类进行分割，可讲整个矩阵划分为：

$K^d=\left(\begin{array}{ccc}{K}_{\mathrm{ii}}^d& K_{\mathrm{if}}^d& K_{\mathrm{ic}}^d\\ K_{\mathrm{fi}}^d& K_{\mathrm{ff}}^d& K_{\mathrm{fc}}^d\\ K_{\mathrm{ci}}^d& K_{\mathrm{cf}}^d& K_{\mathrm{cc}}^d\end{array}\right)$

将Boundary mode乘在stiffness矩阵的左右两边，可得到：

$L^d={\mathrm{\Φ}}^{\mathrm{dT}}{K}^d{\mathrm{\Φ}}^d=\left(\begin{array}{cc}{L}_{\mathrm{ff}}^d& L_{\mathrm{fc}}^d\\ L_{\mathrm{fc}}^{\mathrm{dT}}& L_{\mathrm{cc}}^d\end{array}\right)$

$L_{\mathrm{ff}}^d={\mathrm{\Φ}}_{\mathrm{if}}^{\mathrm{dT}}{K}_{\mathrm{ii}}^d{\mathrm{\Φ}}_{\mathrm{if}}^d+K_{\mathrm{if}}^{\mathrm{dT}}{\mathrm{\Φ}}_{\mathrm{if}}^d+{\mathrm{\Φ}}_{\mathrm{if}}^{\mathrm{dT}}{K}_{\mathrm{if}}^d+K_{\mathrm{ff}}^d$

$L_{\mathrm{fc}}^d={\mathrm{\Φ}}_{\mathrm{if}}^{\mathrm{dT}}{K}_{\mathrm{ii}}^d{\mathrm{\Φ}}_{\mathrm{ic}}^d+K_{\mathrm{if}}^{\mathrm{dT}}{\mathrm{\Φ}}_{\mathrm{ic}}^d+{\mathrm{\Φ}}_{\mathrm{if}}^{\mathrm{dT}}{K}_{\mathrm{ic}}^d+K_{\mathrm{fc}}^d$

$L_{\mathrm{cc}}^d={\mathrm{\Φ}}_{\mathrm{ic}}^{\mathrm{dT}}{K}_{\mathrm{ii}}^d{\mathrm{\Φ}}_{\mathrm{ic}}^d+K_{\mathrm{ic}}^{\mathrm{dT}}{\mathrm{\Φ}}_{\mathrm{ic}}^d+{\mathrm{\Φ}}_{\mathrm{ic}}^{\mathrm{dT}}{K}_{\mathrm{ic}}^d+K_{\mathrm{cc}}^d$

将这些项替换到原FETI-DP方法中，可得到降维后的系统：

$\left(\begin{array}{cc}{G}_{\mathrm{ff}}& G_{\mathrm{fc}}\\ G_{\mathrm{fc}}^T& G_{\mathrm{cc}}\end{array}\right)\left(\begin{array}{c}\mathrm{\λ}\\ q_c\end{array}\right)=\left(\begin{array}{c}{h}_f\\ h_c\end{array}\right)$

$G_{\mathrm{ff}}=\mathrm{\Σ}{B}_f^d{L_{\mathrm{ff}}^d}^{-1}{B}_f^{\mathrm{dT}}$

$G_{\mathrm{fc}}=\mathrm{\Σ}{B}_f^d{L_{\mathrm{ff}}^d}^{-1}{L}_{\mathrm{fc}}^d{B}_c^d$

$G_{\mathrm{cc}}=\mathrm{\Σ}{B}_c^{\mathrm{dT}}(L_{\mathrm{fc}}^{\mathrm{dT}}{L_{\mathrm{ff}}^d}^{-1}{L}_{\mathrm{fc}}^d-L_{\mathrm{cc}}^d)B_c^d$

$h_f=\mathrm{\Σ}{B}_f^d{L_{\mathrm{ff}}^d}^{-1}({\mathrm{\Φ}}_{\mathrm{if}}^{\mathrm{dT}}{f}_i^d+f_f^d)$

$h_c=\mathrm{\Σ}(B_f^{\mathrm{dT}}{L}_{\mathrm{fc}}^{\mathrm{dT}}{L_{\mathrm{ff}}^d}^{-1}({\mathrm{\Φ}}_{\mathrm{if}}^{\mathrm{dT}}{f}_i^d+f_f^d)-B_c^{\mathrm{dT}}({\mathrm{\Φ}}_{\mathrm{ic}}^{\mathrm{dT}}{f}_i-f_c^d))$

求解该系统得到一组向量解u',带入公式u＝Φu'即可得到顶点的位移。

由顶点位移计算应力的公式为：

σ＝DBu

将上面得到的顶点位移量带入该式，即得到网格模型的应力分布。

4.3.2)应力分布情况的可视化：本发明采用Von-Mises应力作为应力分布可视化所使用的值。Von-Mises应力的计算公式为：

${\mathrm{\σ}}_{\mathrm{vm}}=\frac{1}{\sqrt{2}}({({\mathrm{\σ}}_x-{\mathrm{\σ}}_y)}^2+{({\mathrm{\σ}}_y-{\mathrm{\σ}}_z)}^2+{({\mathrm{\σ}}_z-{\mathrm{\σ}}_x)}^2+6({\mathrm{\τ}}_{\mathrm{xy}}^2+{\mathrm{\τ}}_{\mathrm{yz}}^2+{\mathrm{\τ}}_{\mathrm{zx}}^2))$

其中σ_x,σ_y,σ_z,τ_xy,τ_yz,τ_zx为上一步骤中算得应力σ的6个分量。

在算得Von-Mises应力值后，本发明将该值的分布采用颜色插值的方式显示给用户作为应力分布情况的估算(如图3所示)。

5)应力分布情况的准确计算：本发明采用FETI-DP方法(FARHAT,C.,LESOINNE,M.,LETALLEC,666P.,PIERSON,K.,AND RIXEN,D.2001.Feti-dp:a dual-primal unified fetimethod part i:A faster alternative to the two-level feti method.International Journal for NumericalMethods in Engineering50,7,1523–1544)，计算网格模型应力分布的准确值。在使用该方法求得顶点位移的准确值之后，采用和前述步骤一样的方法，将应力分布转换为Von-Mises应力分布，并可视化。

进一步地说，一种为用户提供网格模型形状编辑工具，且在编辑过程中即能快速估算应力分布情况以辅助用户进行编辑操作，又能准确求解应力分布情况的系统，包括以下步骤：

a)从输入四边形网格模型自动或手工抽取骨骼结构；

b)根据骨骼结构将四边形网格划分为区域；

c)用3)中的装置预计算每个域的参数化stiffness矩阵。

d)用2)中的装置预计算Boundary mode估算数据。

e)用4.1)中的装置对网格进行编辑，修改网格形状。

f)用4.2)中的装置更新应力求解系统中的运算数据。

g)用4.3)中的方法求解应力系统，将结果以色彩或形变的方式显示，辅助用户进行编辑。

h)用原FETI-DP方法求解应力分布的准确值，将结果以色彩或形变的方式显示。

实施例：

发明人在一台配备Intel i73770k(4核3.4GHz)，32GB内存的机器上实现了本发明的一个实施实例。发明人所采用的材料物理系数为近似ABS塑料的物理系数，其杨氏模数为2.3e9，密度为1040千克/立方米，泊松系数为0.35。图2中所示的鸟类模型有50995个顶点，26020个四面体元。在网格模型预处理阶段，该模型被划分为21个网格域。以该模型作为输入，为估算实时估算应力分布而建立的求解系统，平均每个网格域约有240个自由度。编辑过程中系统的更新时间平均为1.36秒，求解时间平均为416.8毫秒。为计算准确应力分布情况而建立的求解系统，平均每个网格域有7285个自由度，计算时间为17秒左右。计算所得应力分布情况可见图3。

发明人对比了本方法与不采用区域划分技术而整体求解的方法之间的效率差异。以上文所述鸟类模型作为输入，若不采用区域划分技术，整体对网格模型进行应力分布求解时，所建立的应力求解系统的求解时间约为123.5秒，本发明所采用的方法可带来7倍左右的加速。本发明所开发的估算算法更是可实时地预览应力分布情况。

发明人还将编辑过的模型实际打印出来，验证本发明所算得应力分布的准确性。发明人在所实现的实施实例上对网格模型进行编辑，并记录下输入模型在不同受力情况下的应力分布情况。在实际打印出来的打印结果上，发明人施加同等大小的真实拉力，并观察打印结果的形变情况。实验证明，打印结果可以承受软件模拟中被认为安全的受力情况，分析结果安全可靠。

除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。

Claims (7)

1.一种面向三维制造技术的快速应力分析方法，包括以下步骤：

1)对输入网格模型的预处理；

2)Boundary mode近似参数函数曲面系数的预计算；

3)参数化stiffness矩阵的预计算；

4)在用户进行编辑时即时更新应力求解系统的运算数据并实时估算应力分布情况；

5)编辑完成后提供应力分布情况的准确计算。

2.根据权利要求1所述的面向三维制造技术的快速应力分析方法，其特征是：所述2)包含以下步骤：

2.1)对于每一个网格域，在其s_x,s_y,s_z三个缩放参数上进行采样，形成一组采样点；

2.2)对于每一个网格域，在各个采样点构建stiffness矩阵，然后分别计算采样点处的Boundary mode；

2.3)以s_x,s_y,s_z三个缩放参数作为坐标轴，将Boundary mode矩阵的每一个元素的值视为一个高维曲面，采用最小二乘方法，用一个参数曲面函数逼近这个高维曲面，保存每个元素得到的曲面的系数；

2.4)当需要估算某组缩放参数所对应的Boundary mode时，将这组缩放参数带入到Boundary mode每一项的曲面函数中，即得到Boundary mode的逼近值。

3.根据权利要求1所述的面向三维制造技术的快速应力分析方法，其特征是：所述3)包含以下步骤：

3.1)将stiffness矩阵的每一项分别表示成多项式形式，多项式的参数为网格域的缩放参数以及材料的物理参数；

3.2)遍历四面体网格的拓扑关系，为矩阵的每一项的多项式系数分配内存空间；

3.3)遍历四面体网格的顶点坐标，预计算矩阵每一项的多项式系数；

3.4)当网格因为用户的编辑进行缩放或者其物理参数发生改变时，将编辑后的缩放系数和物理参数当作多项式参数带入3.2)和3.3)中预计算得到的多项式，即可得到stiffness矩阵。

4.根据权利要求1所述的面向三维制造技术的快速应力分析方法，其特征是：所述4)包含以下步骤：

4.1)编辑操作；

4.2)即时更新应力求解系统的运算数据；

4.3)应力分布情况的快速估算。

5.根据权利要求4所述的面向三维制造技术的快速应力分析方法，其特征是：所述4.1)中将剖分得到的网格域分成两类：刚体运动域和可形变域，可形变域的顶点坐标值在每一个所邻接的刚体运动域中保存一份变换到刚体域局部坐标系中的局部坐标值以及对应于该刚体运动域的混合权重，对于刚体运动域，用户可进行平移，旋转以及缩放操作，当可形变网格域所邻接的刚体运动域经编辑后，重新计算可形变域的顶点坐标：将每个邻接的刚体运动域中所保存的形变域的顶点坐标值按照刚体运动域的当前局部坐标系重新变换到可形变域的局部坐标系中，乘以权重并累加。

6.根据权利要求4所述的面向三维制造技术的快速应力分析方法，其特征是：所述4.2)中当用户对某一个刚体运动域进行平移编辑时，不需要对应力求解系统进行更新；当用户对一个刚体运动域进行旋转编辑时，根据编辑后的域局部坐标系重新计算该域在应力求解系统中对应的乘子矩阵；当用户对一个刚体运动域进行缩放操作时，用重新计算该域在应力求解系统中所对应的stiffness矩阵，重新估算该域在应力求解系统中所对应的Boundarymode；当可形变域因为邻接的刚体运动域发生受到编辑影响而受到影响时，重新计算该域的stiffness矩阵。

7.根据权利要求4所述的面向三维制造技术的快速应力分析方法，其特征是：所述4.3)中用Boundary mode对所有刚体运动域的stiffness矩阵进行降维，将解空间的维度降至Boundary mode所能表达的形变空间，对降维后的矩阵进行原FETI-DP算法所进行的运算，求得原问题在Boundary mode形变空间内所能表达的一个逼近解。