CN104794742A - Balloon expansion animation simulation method based on finite element method - Google Patents

Abstract

The invention discloses a balloon expansion animation simulation method based on a finite element method. The method comprises the following steps of: setting intensity of pressure of gas, balloon density, and balloon initial thickness, reading data from an obj model as un-deformed shape of the balloon, wherein the data comprise junction position and triangle information; performing pre-computation on the model, namely computing the quality of each junction, initial volume of the balloon and the like; calculating out each junction elasticity and gas pressure based on the finite element method, and calculating out the position and speed of each junction at next time by Newton-Raphson iteration method; correcting the calculated out the new position and speed by a volume keeping algorithm, in the way, the volume of the balloon is accordant with the initial volume; rendering the model by POV-Ray software as the image of the next frame; performing the calculation of the next frame by using the new result as the current position, until the balloon is stable and does not deform again. The method disclosed by the invention is simple and efficient, is applicable to deformation of a volume-keeping thin-shell model; and the relatively real expansion effect of the balloon and different effects of balloon with different thicknesses during deformation can be realized.

Description

A kind of ballooning animation simulation method based on Finite Element Method

Technical field

The invention belongs to field of Computer Graphics, be specifically related to a kind of ballooning animation simulation method based on Finite Element Method.

Background technology

Along with the high speed development of computing machine, the simulation of physically based deformation is used in the aspects such as game, animation, video display special efficacy in a large number, the distortion of the deformable body such as motion, muscle skin of the fluid such as motion, current of such as rigid body, the dynamic similation of cloth, hair etc.In the simulation of deformable body, the simulation of Thin shell model is a large difficult point, and balloon is typical Thin shell model, and is the Thin shell model that deformation is larger.

Thin shell model and sheet model closely similar, they are the objects that thickness is all far smaller than other two dimensions.The difference of Thin shell model and sheet model is the shape of non-deformation state.Thin shell model is curved surface when not being out of shape, such as leaf, cap, balloon etc.And sheet model is plane when not being out of shape, such as cloth.The thickness dimension of Thin shell model and sheet model is far smaller than other dimensions, so the general three-dimensional body cannot be used as in space is simulated, otherwise easily occurs ill-conditioning problem, greatly reduces the degree of accuracy of simulation.Due to the similarity of Thin shell model and sheet model, a lot of people makes to use the same method and to simulate two kinds of models, ignores thickness, model is approximately plane in space or curved surface calculates.But for Thin shell model, particularly very irregular Thin shell model, be difficult to calculate its two-dimensional parameter coordinate, the approximation method therefore has larger difficulty when realizing.On the other hand, ignore thickness and lost very important information, the balloon model in such as this method for some object.Different thickness, stretching that balloon produces everywhere is different, and display effect also can be different, ignore thickness information, untrue when obviously can cause drafting.

In addition, balloon, compared with other Thin shell models, has unique characteristic: what (1) balloon played a major role in distortion is the power produced that stretches, and the bending power produced is very little comparatively speaking.(2) balloon is in distortion, square section does not have substantially obvious mistake cut, and namely at thickness direction, does not have obvious mistake to cut.(3) balloon is in deformation process, keeps constancy of volume.Here volume refers to the volume of balloon itself, and does not comprise the volume of gas in balloon.Along with the expansion of balloon, the lower thickness of balloon, vice versa.The characteristic of these uniquenesses further increases the difficulty of balloon deformation simulative, and current method cannot realize these characteristics.

Summary of the invention

The object of the present invention is to provide a kind of ballooning animation simulation method based on Finite Element Method, be applied in ballooning simulation animation, simulate effect is true, efficient, and can realize the different distortion effects under different-thickness, be that prior art cannot realize.

For achieving the above object, the present invention adopts following technical scheme:

Based on a ballooning animation simulation method for Finite Element Method, the steps include:

(1) gas pressure intensity of balloon interior, balloon density, balloon original depth parameter are set, and from obj model file, read in grid data as the undeformed shape of balloon, grid data comprises a series of triangle, and each leg-of-mutton data comprise leg-of-mutton three site positions and triangle information representing;

(2) precomputation is carried out to the model read in step (1), comprise the quality, the initial volume of balloon, the initial volume of gas that calculate each node;

(3) based on Finite Element Method, the result in conjunction with precomputation calculates elastic force, the gaseous tension of each elementary cell, thus obtains each node elastic force, gaseous tension, and uses Newton-Raphson alternative manner to obtain each node in the position of subsequent time and speed;

(4) application protects volume algorithm correction step (3) each site position of obtaining and speed as each node in the position of next frame and speed, and the volume of balloon and initial volume are consistent;

(5) POV-Ray software is used to play up, as the image of next frame the position of each node at next frame;

(6) using the position of each node obtained in step (4) and speed as current location and speed, get back to the calculating that step (3) carries out the position of each node of next frame, until balloon settles out not redeformation.

In described step (1), leg-of-mutton area and distribution are as far as possible evenly; The model read in from obj file, as the undeformed shape of balloon, suitably carries out convergent-divergent to obj model, makes subsequent calculations degree of accuracy higher.

Described step (3) elementary cell is tri-prism element, tri-prism element be triangular element conventional in additive method basis on increase thickness information, namely each tri-prism element comprises 3 nodes, 1 thickness scalar; The original depth of each tri-prism element is the original depth of balloon, in calculating afterwards, Finite Element Method upgrades the position of 3 nodes, protects in volume algorithm and upgrades the position of node and the thickness of each tri-prism element, ensure that the volume of balloon is constant in deformation process with this.

In described step (2), the quality of each node calculates according to each tri-prism element volume and density and mass distribution weight; After obtaining the quality of each tri-prism element, discrete on 3 nodes according to weight.

In described step (2), the initial volume of balloon is obtained by each tri-prism element volume integration.

In described step (2), the tetrahedron integration that the initial volume of gas is formed by tri-prism element base triangle and initial point obtains.

In described step (2), also comprise matrix when calculating non-deformation state in deformation gradient in precomputation, calculate in advance and store, facilitate follow-up use.

In described step (3), obtain each balloon node elastic force, gaseous tension based on Finite Element Method, and use Newton-Raphson alternative manner to obtain each balloon node in the position of subsequent time and speed flow process to be:

(31) in each tri-prism element, the elastic force that calculating crunode is subject to and gaseous tension, in the corresponding elastic force vector sum gaseous tension vector of the node that is added to;

(32) the node elastic force vector sum pressure vector tectonic system linear equation utilizing step (31) to calculate;

(33) adopt method of conjugate gradient to solve linear equation, in solution procedure, obtain the variation of elastic force and pressure as required, as stiffness matrix and matrix of coefficients and vectorial product;

(34) result obtained according to step (33) upgrades site position and speed, gets back to step (31) and enters next iteration; Iterations is 3-5 time, determines according to model; The result of last iteration is the net result of step in claim 1 (3).

The calculating of the elastic force vector in described step (32) obtains site position according to current node position and when not being out of shape, its calculation process is: to each triangular prism, calculate deformation gradient, strain tensor, stress tensor, deformation gradient to the derivative of 3 nodes, then to calculate in each tri-prism element elastic force that 3 nodes are subject to and be added in the elastic force vector of node; Finally, the elastic force vector of node is total elastic force that this node is subject to.

In described step (32), the calculating of gaseous tension vector obtains according to current gas volume, initial gas volume, initial gas pressure, tri-prism element base area and normal vector, calculation process is: the volume first calculating current gas, then current gas pressure is obtained according to current gas volume and initial gas volume and pressure, then obtain the gaseous tension that each tri-prism element is subject to, be finally added in the pressure vector of node by discrete for tri-prism element gaseous tension to 3 nodes; Finally, the pressure vector of node is total gaseous tension that this node is subject to.

The computing method of the variation of described step (33) elastic force are: obtain the variation of deformation gradient, the variation of strain tensor, the variation of stress tensor, then calculate the elastic force variation of 3 nodes in each tri-prism element and are added in the elastic force variation vector of node.

Each node (whose) position that in described step (4), application guarantor's volume algorithm correction step (3) is obtained and speed flow process are:

(41) thickness information of each tri-prism element is converted into virtual nodes, as the 4th node of tri-prism element, and by the node list of this node join system;

(42) in each tri-prism element, the new quality adding totally 4 nodes of the 4th node is calculated;

(43) calculate the objective function protecting volume, i.e. the difference of current balloon volume and initial volume, then calculating target function comprises the derivative of virtual nodes to each node;

(44) calculate the position correction amount of each node according to the above results, comprise virtual nodes;

(45) according to correction, upgrade the position and the speed that comprise all nodes of virtual nodes, then in each tri-prism element, calculate new thickness according to the reposition of 4 nodes.

The advantage that the present invention compared with prior art has and good effect as follows:

(1) the present invention is based on Finite Element Method, propose the deformation algorithm of simple Thin shell model, compared with prior art simple efficient, and don't can thickness information be lost.

(2) the present invention is protecting volume algorithm based on the deformation algorithm basis of finite element adds, and realize the Thin shell model modeling algorithm protecting volume, prior art does not realize real guarantor's volume effect.

(3) the present invention is by algorithm application in ballooning simulation animation, and simulate effect is true, efficient, and can realize the different distortion effects under different-thickness, be that prior art cannot realize.

Accompanying drawing explanation

Fig. 1 is based on the ballooning animation simulation method flow diagram of Finite Element Method;

Fig. 2 finite element method process flow diagram;

Fig. 3 protects volume algorithm flow chart.

Embodiment

In order to make the those skilled in the art of the art better understand the present invention, below in conjunction with accompanying drawing, a kind of ballooning animation simulation method based on Finite Element Method provided by the present invention being described in further detail, but not being construed as limiting the invention.

The present invention is using simple tri-prism element as the elementary cell of Finite Element Method, thickness information is introduced, the distortion of this rubber Thin shell model of simulation balloon, realize the vital role of thickness in balloon distortion, and carry out analog gas with simple uniform gas model, realize the dynamic process blown up a balloon, the different distortion process of balloon when inflating of simulation different-thickness.The method speed is fast, and substantially realize the characteristic of balloon, the volume of balloon substantially remains unchanged in deformation process, and deformation process meets the deformation effect of true balloon substantially, and the distortion of the balloon of different-thickness is different, is that additive method institute is inaccessiable.

Fig. 1 is the ballooning animation simulation method flow diagram based on Finite Element Method.The key step of method is as follows:

(1) parameters, as gas pressure intensity, balloon density etc.

(2) from obj model, read data, require that model is triangular mesh, as far as possible evenly.

(3) precomputation.Comprise calculating crunode quality, balloon initial volume etc.

(4) position and the speed of each node of finite element method next frame is used.

(5) site position using guarantor's volume algorithm correction previous step to obtain and speed, make balloon volume equal with initial volume.

(6) POV-Ray software rending model is used, as the image of next frame.

(7) get back to step (4), using the result that calculates as current location, continue the image calculating next frame, until balloon settles out not redeformation.

Step (4) is of the present invention one large pith, finite element method part.This part is a longer part consuming time in algorithm, and the speed of this part determines the arithmetic speed of total algorithm, and the rough location of the node obtained determines simulate effect of the present invention.As shown in Figure 2, method is as follows for concrete solution procedure:

(1) in each cell, the elastic force that calculating crunode is subject to and gaseous tension, in the corresponding vector of the node that is added to.

(2) the node elastic force vector sum pressure vector tectonic system linear equation using previous step to calculate.

(3) use method of conjugate gradient to solve linear equation, in solution procedure, obtain the variation of elastic force and pressure as required, as stiffness matrix and matrix of coefficients and vectorial product.

(4) result obtained according to previous step upgrades site position and speed, gets back to step (1) and enters next iteration.Iterations is generally 3-5 time, determines according to model.The result of last iteration is the net result of step (4) in Fig. 1 process flow diagram.

In Fig. 1 process flow diagram, step (5) protects volume algorithm is another large major part of the present invention, and this part revises the position of node further on the basis of Finite Element Method, makes it meet the characteristic of balloon constancy of volume.This part determines the effect whether the present invention can simulate different-thickness.As shown in Figure 3, method is as follows for concrete solution procedure:

(1) thickness information of each unit is converted into virtual node, and this node is added in the node list of system.

(2) the new quality of 4 nodes in each cell, is calculated.

(3) objective function protecting volume is calculated, i.e. the difference of current balloon volume and initial volume.Then calculating target function is to the derivative of each node, comprises virtual nodes.

(4) calculate the position correction amount of each node according to the above results, comprise virtual nodes.

(5) according to correction, the position and the speed that comprise all nodes of virtual nodes is upgraded.Then in each cell, new thickness is calculated according to the reposition of 4 nodes.

Lower mask body describes in detail, the steps include:

1) model importing, optimum configurations and precomputation

Require that model uses triangular mesh, model data file is .obj form, and Triangle-Profile as far as possible evenly.Read in the file of .obj form, read the information of node and leg-of-mutton information, node initial position is set namely put the coordinate of i at the position in 0 moment i.e. file mid point i, initial velocity is set be 0, parameters initial gas pressure p ₀balloon density p and initial balloon thickness h.Consider the characteristic of balloon and Thin shell model, the triangular prism that each unit is made up of leg-of-mutton three nodes and thickness profile, while introducing thickness information, and don't the complexity of excessive increase method.Volume and density calculation node quality according to unit:

$m_i=\underset{{e\∈T}^i}{\mathrm{\Σ}}{w}_i\mathrm{\ρ}\·V^e$

Wherein Τ ⁱrepresent the unit set comprising node i, V ^ebeing the volume of tri-prism element, is the product of floorage and thickness.W _ibe the weight that quality is loaded on node i by unit e, use 1/3 here, in follow-up calculating, the quality that different weight calculation is new can be used.

Calculate initial gas volume, the volume of gas is the volume of balloon interior, can be obtained by triangular mesh, and the initial point of each triangle and coordinate system can form a tetrahedron, all tetrahedral volumes are added, are the volume of balloon interior i.e. the volume of gas:

$V_0=\underset{e}{\mathrm{\Σ}}\frac{1}{3}\·\frac{1}{3}({\stackrel{\‾}{x}}_1^e+{\stackrel{\‾}{x}}_2^e+{\stackrel{\‾}{x}}_3^e)\·n^e\·A^e$

Wherein, A ^ethe volume of unit base triangle, n ^ethe normal vector of unit base triangle, for three nodes of tri-prism element.The center of three nodes is tetrahedral height in the projection in normal vector direction, is tetrahedral volume with 1/3rd of the product of the area of base triangle.

For convenience of subsequent calculations, this step does following precomputation to each unit:

$M^e={\left[\begin{array}{ccc}{\stackrel{\‾}{e}}_{12}& {\stackrel{\‾}{e}}_{13}& \stackrel{\‾}{d}{\stackrel{\‾}{h}}^e\end{array}\right]}^{-1}$

$\stackrel{\‾}{d}=\frac{{\stackrel{\‾}{e}}_{12}\×{\stackrel{\‾}{e}}_{13}}{|{\stackrel{\‾}{e}}_{12}\×{\stackrel{\‾}{e}}_{13}|}$

Wherein in unit, 3 nodes of base triangle are thickness is limit for the unit normal vector of unit base triangle.

2) elastic force that is subject to of calculating crunode

From this step, enter iteration, each iteration all calculates speed and the position of each node of next moment according to the speed of current node and position, i.e. the state of next frame in animation.When each iteration starts, the node speed of the subsequent time obtained in use last iteration and position are as present speed and position.If first time iteration, then use step 1) in node initial position and initial velocity as current location and speed.

Calculate the deformation gradient of each unit.If the current location of 3 nodes is in unit e element thickness is h ^e∈ R, limit d is the unit normal vector of unit base triangle, then by the deformation gradient of following formula computing unit e:

F ^e＝[e ₁₂e ₁₃dh ^e]·M ^e

$d=\frac{e_{12}\×e_{13}}{|e_{12}\×e_{13}|}$

Wherein, M ^estep 1) in precomputation result, represent the information of non-deformation state unit.Deformation gradient represents the deformation degree of object, relevant with non-deformation position with the position of current node, and current location is larger with non-deformation position gap, then illustrate that the deformation of object is larger.The form of this deformation gradient is derived from tri-prism element, according to any point in triangular prism before being deformed be out of shape after coordinate obtain.Different unit has different deformation gradient expression-forms, and the complicated and simple of this form plays an important role to follow-up calculating.

Obtain the strain tensor of unit according to deformation gradient, strain tensor represents the distortion of object, more can reflect the characteristic of object, and this method uses Green tensor:

$E^e=\frac{1}{2}(F^{\mathrm{eT}}{F}^e-I)$

The 1st Piola-kirchhoff stress tensor of each unit is calculated according to strain tensor.After object generation deformation, interior of articles can produce stress, as producing the elastic force of opposing stretching after spring elongation.Stress tensor is the tensor form of stress.Different materials, the form of stress tensor is different, with the relating to parameters of properties of materials and material.This method is common to different materials, but the singularity of large deformation due to the elastomeric material of balloon, this method is mainly for nonlinear material.The strain tensor of St.Venant-Kirchhoff material is calculated by following formula:

P ^e(F ^e)＝F ^e[2μE ^e+λtr(E ^e)I ^e]

Wherein, μ, λ are the Lame Coefficient of material.

Finally, calculate according to stress tensor the elastic force that each node is subject to:

$f_i^{\mathrm{elastic}}=\underset{e}{\mathrm{\Σ}}{f}_i^e=\underset{e}{\mathrm{\Σ}}-V^{e}P\left(F^e\right)\frac{{\∂F}^e}{{\∂x}_i}$

In this step, first according to the deformation gradient of unit e obtained before to the derivative of three nodes, then obtain according to above formula the elastic force that in this unit, three nodes are subject to, be added in the elastic force vector of 3 nodes.All above-mentioned calculating is carried out to all unit, finally can obtain the elastic force that each node is subject to.Must, the elastic force that node is subject to is exactly the summation of the elastic force that it is subject in all unit around.

Sum up the computation process of this step: to each unit, calculate deformation gradient, strain tensor, stress tensor, deformation gradient to the derivative of 3 nodes, then to calculate in each unit elastic force that 3 nodes are subject to and be added in the elastic force vector of node.Finally, the elastic force vector of node is total elastic force that this node is subject to.

3) gaseous tension that is subject to of calculating crunode

In the gas replenishment process of balloon, the flowing of gas is not obvious, and gas and balloon reach equilibrium state very soon, and therefore this method adopts simple uniform gas model, and namely in balloon, the pressure of each point is identical.

In the balloon closed, the pressure of gas and the product of volume are constants.First, with step 1) the same, calculate the current volume of gas:

$V=\underset{e}{\mathrm{\Σ}}\frac{1}{3}\·\frac{1}{3}(x_1^e+x_2^e+x_3^e)\·n^e\·A^e$

The current pressure of gas can be obtained by current volume, initial volume, initial pressure:

$p=\frac{p_0{V}_0}{V}$

Pressure of each point of gas interior is all equal, and therefore according to the area of pressure and unit base triangle, can obtain the size of gaseous tension suffered by unit, the direction of gaseous tension is the direction of the normal vector of base triangle.The gaseous tension that node is subject in this unit is 1/3rd of the gaseous tension that unit is subject to.Calculate the gaseous tension of all unit, and it discretely finally can be obtained the gaseous tension of node on unit node, namely the gaseous tension of node is the superposition of gaseous tension on node that around it, all unit are subject to.

$f_i^p=\underset{{e\∈T}^i}{\mathrm{\Σ}}\frac{1}{3}{\mathrm{pA}}^e{n}^e=\underset{{e\∈T}^i}{\mathrm{\Σ}}\frac{1}{3}p\·\frac{1}{2}(e_{12}^e\×e_{13}^e)$

4) tectonic system linear equation solving

If first time iteration, then establish initial value be defined as follows variable:

${\mathrm{\Δx}}_{\left(k\right)}=x_{(k+1)}^{n+1}-x_{\left(k\right)}^{n+1},{\mathrm{\Δv}}_{\left(k\right)}=v_{(k+1)}^{n+1}-v_{\left(k\right)}^{n+1}$

Then basis formula construction system linear equation below, uses Newton-Raphson method iterative, obtains a series of approximate value and as k → ∞, $\underset{k\→\∞}{\lim }{x}_{\left(k\right)}^{n+1}=x^{n+1},\underset{k\→\∞}{\lim }{v}_{\left(k\right)}^{n+1}=v^{n+1}.$ This process repeats can obtain result more accurately 3-5 time, determine according to concrete model, in each iteration, use conjugate gradient method to obtain the solution of linear equation, after obtaining solution, get back to step 2) enter next iteration, enter step 5 after iteration terminates).

A△v _(k)＝b

$A=M+\mathrm{\ΔtC}+{\mathrm{\Δt}}^2(K^{\mathrm{elastic}}\left(x_{\left(k\right)}^{n+1}\right)+K^p\left(x_{\left(k\right)}^{n+1}\right))$

$b=M(v^n-v_{\left(k\right)}^{n+1})+\mathrm{\Δt}(f^{\mathrm{elastic}}\left(x_{\left(k\right)}^{n+1}\right)+f^{\mathrm{ext}}+f^p\left(x_{\left(k\right)}^{n+1}\right)-{\mathrm{Cv}}_{\left(k\right)}^{n+1})$

$v_{(k+1)}^{n+1}={\mathrm{\Δv}}_{\left(k\right)}+v_{\left(k\right)}^{n+1}$

$x_{(k+1)}^{n+1}=\mathrm{\Δt}\·{\mathrm{\Δv}}_{\left(k\right)}+x_{\left(k\right)}^{n+1}$

Wherein, M is mass matrix, is diagonal matrix, and diagonal line is the quality of each node, and other elements are 0.△ t is time step, i.e. the mistiming of current time and subsequent time.V ⁿbe the rate matrices of moment n, C is damping matrix, generally uses Rayleigh damping, C=α M+ β K ^elastic, f ^extother power except elastic force, pressure, such as gravity.K ^elasticand K ^pbe the stiffness matrix of elastic force and the matrix of coefficients of gaseous tension respectively, they are power derivatives to node coordinate, that is:

$K^{\mathrm{elastic}}=-\frac{{\∂f}^{\mathrm{elastic}}}{\∂x},K^p=-\frac{{\∂f}^p}{\∂x}$

In method before, mostly directly obtain this two matrixes.In this method, in order to accelerate arithmetic speed and reduce the complexity of computing, and consider that these two matrixes are not occur separately in solution procedure, but occur with the form of vector product, i.e. K ^elasticw, K ^pthe form of w occurs, w is any vector, and therefore our demand goes out the value of product.

For the stiffness matrix of elastic force, to any vectorial w, make δ x=-w, then solving target is δ f ^elastic=-K ^elasticδ x

Concrete solution procedure is as follows:

To each unit:

● obtain the variation δ F of deformation gradient ^e=[δ e ₁₂δ e ₁₃δ dh ^e] M ^e

● obtain corresponding strain tensor

● obtain corresponding stress tensor δ P (F; δ F)=δ F [2 μ E+ λ tr (E) I]+F [2 μ δ E+ λ tr (δ E) I]

● obtain the elastic force differential of each node in unit and in the elastic force differential vector of the node that is added to.

After all unit solve end, required by the elastic force differential vector of node is.

The matrix of coefficients solution procedure of gaseous tension is as follows:

● to each node i, obtain the normal vector of node

● by the normal vector of all nodes and corresponding δ x _i=-w _icalculate inner product and be added, obtaining scalar

● to each node i, obtain the Pressure differential of node the concrete derivation of this formula does not repeat here.

It is more than the solution procedure of elastic force differential and Pressure differential.In the solution procedure of system linear equation, as long as run into K ^elasticw, K ^pthe form of w just solves by above-mentioned steps.

5) volume algorithm is protected

Previous step obtains node behind the position and speed in next moment, may cause the volume change of balloon, and therefore this step uses the position and speed of protecting volume algorithm correction node, makes the constancy of volume of balloon.

In order to shortcut calculation, first the amount of thickness h of each unit is converted into virtual point, as the 4th node of unit.This node is positioned at the center of triangular prism upper bottom surface, and other three nodes are at the bottom surface of triangular prism.

$x_4^e=\frac{1}{3}(x_1^e+x_2^e+x_3^e)+{\mathrm{dh}}^e$

Recalculate the quality of each node according to new node, original three node Mass Calculation weights are 1/6, and the node quality weight newly increased is 1/2.Then the objective function protecting volume constraint is calculated, i.e. the difference of current balloon volume and initial volume the calculating of volume is still obtained by the volume integral of each tri-prism element.Then calculating target function is to the derivative of each node:

$\frac{\∂C\left(x\right)}{{\∂x}_i}=\underset{e}{\mathrm{\Σ}}\frac{{\∂C}^e\left(x\right)}{{\∂x}_i}=\underset{e:t_1^e=i}{\mathrm{\Σ}}\frac{1}{2}{e}_{24}\×e_{23}+\underset{e:t_2^e=i}{\mathrm{\Σ}}\frac{1}{2}{e}_{13}\×e_{14}+\underset{e:t_3^e=i}{\mathrm{\Σ}}\frac{1}{2}{e}_{14}\×e_{12}+\underset{e:t_4^e=i}{\mathrm{\Σ}}\frac{1}{2}{e}_{12}\×e_{13}$

This derivative is obtained by the integration of each triangular prism volume to node, wherein it is the index of 4 nodes in unit e.Finally, the correction of each site position is obtained according to above result:

${\mathrm{\Δx}}_i=-\frac{\left({N+N}_1\right)w_i}{\underset{j}{\mathrm{\Σ}}{w}_j}\frac{w_{i}C\left(x\right)}{\underset{j}{\mathrm{\Σ}}{\left|\right|\frac{\∂C\left(x\right)}{{\∂x}_j}\left|\right|}^2}\·\frac{\∂C\left(x\right)}{{\∂x}_i}$

Wherein w _ithe weight of node i in volume constraint, the inverse of service property (quality), i.e. w here _i=1/m _i'.N is original nodal point number, N ₁it is the nodal point number that thickness h converts.Finally upgrade site position according to above correction $x_i^{n+1}=x_i^{n+1}+{\mathrm{\Δx}}_i,$ Speed $v_i^{n+1}=(x_i^{n+1}-x_i^n)/\mathrm{\Δt},$ And virtual nodes is transformed back thickness h.

6) play up

After obtaining site position, use POV-Ray software rending model, display result, then gets back to step 2) carry out the calculating of next frame.

Above a kind of ballooning animation simulation method based on Finite Element Method of the present invention has been described in detail, but obvious specific implementation form of the present invention is not limited thereto.For the those skilled in the art of the art, the various apparent change carried out it when not deviating from spirit and the right of the method for the invention is all within protection scope of the present invention.

Claims (10)

1., based on a ballooning animation simulation method for Finite Element Method, it is characterized in that step is:

(1) gas pressure intensity of balloon interior, balloon density, balloon original depth parameter are set, and from obj model file, read in grid data as the undeformed shape of balloon, grid data comprises a series of triangle, and each leg-of-mutton data comprise leg-of-mutton three site positions and triangle information representing;

(2) precomputation is carried out to the model read in step (1), comprise the quality, the initial volume of balloon, the initial volume of gas that calculate each node;

(3) based on Finite Element Method, the result in conjunction with precomputation calculates elastic force, the gaseous tension of each elementary cell, thus obtains each node elastic force, gaseous tension, and uses Newton-Raphson alternative manner to obtain each node in the position of subsequent time and speed;

(4) application protects volume algorithm correction step (3) each site position of obtaining and speed as each node in the position of next frame and speed, and the volume of balloon and initial volume are consistent;

(5) POV-Ray software is used to play up, as the image of next frame the position of each node at next frame;

(6) using the position of each node obtained in step (4) and speed as current location and speed, get back to the calculating that step (3) carries out the position of each node of next frame, until balloon settles out not redeformation.

2. the ballooning animation simulation method based on Finite Element Method according to claim 1, is characterized in that: in described step (1), and leg-of-mutton area and distribution are as far as possible evenly; The model read in from obj file, as the undeformed shape of balloon, suitably carries out convergent-divergent to obj model, makes subsequent calculations degree of accuracy higher.

3. the ballooning animation simulation method based on Finite Element Method according to claim 1, it is characterized in that: step (3) elementary cell is tri-prism element, tri-prism element be triangular element conventional in additive method basis on increase thickness information, namely each tri-prism element comprises 3 nodes, 1 thickness scalar; The original depth of each tri-prism element is the original depth of balloon, in calculating afterwards, Finite Element Method upgrades the position of 3 nodes, protects in volume algorithm and upgrades the position of node and the thickness of each tri-prism element, ensure that the volume of balloon is constant in deformation process with this.

4. the ballooning animation simulation method based on Finite Element Method according to claim 1, is characterized in that: in described step (2), and the quality of each node calculates according to each tri-prism element volume and density and mass distribution weight; After obtaining the quality of each tri-prism element, discrete on 3 nodes according to weight; In described step (2), the initial volume of balloon is obtained by each tri-prism element volume integration; In described step (2), the tetrahedron integration that the initial volume of gas is formed by tri-prism element base triangle and initial point obtains.

5. the ballooning animation simulation method based on Finite Element Method according to claim 1, it is characterized in that: in described step (2), also comprise matrix when calculating non-deformation state in deformation gradient in precomputation, calculate in advance and store, facilitate follow-up use.

6. the ballooning animation simulation method based on Finite Element Method according to claim 1, it is characterized in that: in described step (3), obtain each balloon node elastic force, gaseous tension based on Finite Element Method, and use Newton-Raphson alternative manner to obtain each balloon node in the position of subsequent time and speed flow process to be:

(31) in each tri-prism element, the elastic force that calculating crunode is subject to and gaseous tension, in the corresponding elastic force vector sum gaseous tension vector of the node that is added to;

(32) the node elastic force vector sum pressure vector tectonic system linear equation utilizing step (31) to calculate;

(33) adopt method of conjugate gradient to solve linear equation, in solution procedure, obtain the variation of elastic force and pressure as required, as stiffness matrix and matrix of coefficients and vectorial product;

(34) result obtained according to step (33) upgrades site position and speed, gets back to step (31) and enters next iteration; Iterations is 3-5 time, determines according to model; The result of last iteration is the net result of step in claim 1 (3).

7. the ballooning animation simulation method based on Finite Element Method according to claim 6, it is characterized in that: the calculating of the elastic force vector in described step (32) according to current node position and be not out of shape time site position obtain, its calculation process is: to each triangular prism, calculate deformation gradient, strain tensor, stress tensor, deformation gradient to the derivative of 3 nodes, then to calculate in each tri-prism element elastic force that 3 nodes are subject to and be added in the elastic force vector of node; Finally, the elastic force vector of node is total elastic force that this node is subject to.

8. the ballooning animation simulation method based on Finite Element Method according to claim 6, it is characterized in that: in described step (32), the calculating of gaseous tension vector is according to current gas volume, initial gas volume, initial gas pressure, tri-prism element base area and normal vector obtain, calculation process is: the volume first calculating current gas, then current gas pressure is obtained according to current gas volume and initial gas volume and pressure, then the gaseous tension that each tri-prism element is subject to is obtained, finally be added in the pressure vector of node by discrete for tri-prism element gaseous tension to 3 nodes, finally, the pressure vector of node is total gaseous tension that this node is subject to.

9. the ballooning animation simulation method based on Finite Element Method according to claim 1, it is characterized in that: the computing method of the variation of described step (33) elastic force are: obtain the variation of deformation gradient, the variation of strain tensor, the variation of stress tensor, then calculate the elastic force variation of 3 nodes in each tri-prism element and be added in the elastic force variation vector of node.

10. the ballooning animation simulation method based on Finite Element Method according to claim 1, is characterized in that:

Each site position that in described step (4), application guarantor's volume algorithm correction step (3) is obtained and speed flow process are:

(41) thickness information of each tri-prism element is converted into virtual nodes, as the 4th node of tri-prism element, and by the node list of this node join system;

(42) in each tri-prism element, the new quality adding totally 4 nodes of the 4th node is calculated;

(43) calculate the objective function protecting volume, i.e. the difference of current balloon volume and initial volume, then calculating target function comprises the derivative of virtual nodes to each node;

(44) calculate the position correction amount of each node according to the above results, comprise virtual nodes;

(45) according to correction, upgrade the position and the speed that comprise all nodes of virtual nodes, then in each tri-prism element, calculate new thickness according to the reposition of 4 nodes.