CN102663816B - Plant leaf wilting simulation method based on physical model - Google Patents

Abstract

The invention discloses a plant leaf wilting simulation method based on a physical model, relating to the field of computer simulation technology. The method comprises the following steps of: constructing a voxel set of a three-dimensional plant leaf model, wherein the voxel set comprises a plurality of hexahedron units used for simulating plant leaf mesophyll internal cells; converting the voxel set into a mass-spring model; after applying an external force to the mass-spring model, simulating the displacement of each peak in the hexahedron units; and according to the displacement of each peak in the hexahedron units, interpolating a coordinate of a peak of a triangular mesh in the plant leaf model. According to the method, the voxel set is constructed through the three-dimensional model, on the basis that the mass-spring model constructed by the voxel set is multi-layered, the model is more close to an internal structure of a real plant leaf, and the connection between mass springs is more stable; an implicit Newmark integration algorithm is employed, simulation can be carried out under a large step length, and the model is ensured to carry out real-time visual simulation.

Description

A kind of plant leaf blade wilting analogy method based on physical model

Technical field

The present invention relates to computer modeling technique field, more particularly to a kind of plant leaf blade wilting mould based on physical model Plan method.

Background technology

The morphosiss modeling and its behavioral trait simulation of plant is that computer graphicss, digital the important of plant field are ground Study carefully content.It is well known that the form of plant is the process of an active development, in varying environment condition and the shadow of anthropic factor Ring, the topological structure of plant can produce significant change., used as the vitals of plant, its form also can be with external rings for leaf Border and its internal physiological process and there is very significant change, such as the factor such as hydropenia, anoxia, high temperature, disease all may draw Play plant leaf blade and the significantly deformation such as curling, wilting occurs.In digital plant field, it is for the intelligent meter of plant characteristic reaches More basic work, it is played to improving numeral plant technology system, the digital plant of enhancing and contacting for agricultural research application Important function.It is contemplated that the wilting analogy method of exploitation plant leaf blade, has great importance.

In recent years, the method for plant leaf blade wilting simulation had obtained certain development, is broadly divided into is become based on skeleton at present Shape and based on the class method of physical deformation two.The method based on framework deformation that Lu et al. is proposed constructs first the leaf vein skeleton of blade (S.Lu, C.Zhao, X.Guo, et al.Venation Skeleton-Based Modeling Plant Leaf Wilting, International Journal of Computer Games Technology, vol.2009, Article ID 890917), skeleton is seen the skeletal chain connected by line segment as, then by setting the moving method of skeleton come driving blade Deformation.This kind of method can obtain more real effect, but the model does not possess any physical significance, completely by user's Interactive operation produces the plant leaf blade form of wilting shape.Chinese patent " the leaves of plants based on physics of Publication No. 101251932 Sub- analogy method " discloses the method based on physics and constructs the physical arrangement of blade first (due to the material feature of blade, employing Mass-spring Model), then realize the dynamic of blade by arranging the methods such as physical attribute, applying external force or change initial condition State deforms, and this method can more reflect the physical attribute of blade than the former, as a result more genuine and believable.But, the method is adopted Mass spring method of attachment is not stable method of attachment；Meanwhile, the integration method that it is adopted all is explicit method, simulation In generally require to arrange very little step-length, cause whole simulation process to need the time (more than 1 hour) grown very much, it is difficult to full Foot carries out the needs of high accuracy or the simulation of large deformation real-time visual.

The content of the invention

(1) technical problem to be solved

The technical problem to be solved in the present invention is：How a kind of plant leaf blade wilting simulation side based on physical model is provided Method, to improve leaf model in mass spring connection stability so that blade wilt simulation effect have more sense of reality.

(2) technical scheme

To solve above-mentioned technical problem, the present invention provides a kind of plant leaf blade wilting analogy method based on physical model, It includes step：

A：The voxel collection of three-dimensional plant leaf model is built, the voxel collection includes multiple for simulating plant blade mesophyll The hexahedral element of internal cell；

B：The voxel collection is converted into into Mass-spring Model；

C：The Mass-spring Model is applied after external force, simulation calculates the position on each summit in the hexahedral element Move；

D：According to the displacement on each summit in the hexahedral element, interpolation goes out the top of the leaf model intermediate cam grid The coordinate of point.

Preferably, step A specifically includes step：

A1：The segmentation of Octree space is carried out to the triangle gridding in leaf model and obtains space Octree；

A2：The spatial surrounding box of the leaf model is calculated, the spatial surrounding box is divided into into multiple hexahedral elements；

A3：All hexahedral elements are traveled through, judges each described hexahedral element with the space Octree successively Subtree space whether intersect, if it is, retaining the current hexahedral element；Otherwise, the current hexahedron is given up Unit；

A4：The traversal institute hexahedral element with a grain of salt, judges successively each described hexahedral element and the subtree Whether the triangle gridding in space intersects, if it is, choosing the current hexahedral element；Otherwise, give up current described Hexahedral element；

A5：The voxel collection of the leaf model is built using all hexahedral elements chosen, and records all The index information on the summit of the triangle gridding included in the hexahedral element chosen.

Preferably, the Mass-spring Model is made up of multiple springs；The spring by the hexahedral element in appoint Two summits of meaning connect and compose, and the summit of the both ends of the spring is denoted as particle.

Preferably, the internal force formula of the spring is as follows：

E=k_s(|x_ij|-L)²；

Wherein, E represents the energy between particle i and particle j, x_ijThe vector difference of particle i and particle j is represented, L represents particle Initial distance between i and particle j, k_sRepresent the stiffness factor of the spring comprising particle i and particle j.

Preferably, the physical kinetics equation of the Mass-spring Model is as follows：

Ma+D (u, v)+R (u)=F；

Wherein, u ∈ R³ⁿRepresent the displacement of particle；v∈R³ⁿFor the first derivative of u, the speed of particle is represented；a∈R³ⁿFor u Second dervative, represent particle acceleration；M∈R^3n×3nRepresent the mass matrix of particle；D (u, v) ∈ R³ⁿRepresent that particle is subject to Damping force；R(u)∈R³ⁿRepresent the internal force that particle is subject to；R represents real number space；N represents particle in the Mass-spring Model Number；F∈R³ⁿRepresent the external force that particle is subject to.

Preferably, the computing formula of internal force R (u) that the particle is subject to is as follows：

Wherein, x_iRepresent the three dimensional space coordinate of particle i.

Preferably, the computing formula of damping force D (u, v) that the particle is subject to is as follows：

D (u, v)=(α M+ ω K (u)) v；

Wherein, α and ω represent the damping parameter of the damping characteristic for controlling the particle spring system；K(u)∈R^3n×3n Tangent stability matrix is represented, computing formula is as follows：

Wherein, x_jThe three dimensional space coordinate of particle j is represented, I represents unit matrix.

Preferably, step C specifically includes step：

C1：Initial setting up is carried out to the stiffness factor of each spring in the Mass-spring Model；

C2：The Mass-spring Model is applied after external force, using implicit expression Newmark integration method to the physical kinetics Equation carries out numerical integration, is calculated the displacement on each summit in the hexahedral element.

Preferably, the implicit expression Newmark integration method is specific as follows：

Input includes：Displacement components u of the particle at the i moment_i, speed v_iAnd acceleration a_i, the i+1 moment be subject to it is outer Power f_i+1, iterationses j_max, and iteration error upper limit TOL, time step t；

Output includes：In the displacement components u at i+1 moment_i+1, speed v_i+1And acceleration a_i+1；

And the implicit expression Newmark integration method includes step：

S1：Setting u_i+1Initial value be u_i, the initial value of variable j is 0；

S2：Whether j judges j less than j from after Jia 1_max, if it is, execution step S3, otherwise, exports current u_i+1、 v_i+1And a_i+1；

S3：According to current u_i+1, calculate the internal force R (u that particle is subject to_i+1), and tangent stability matrix K (u_i+1)；

S4：Calculate Local Damping Matrix C=(α M+ ω K (u))；

S5：Calculate symmetrical matrix A=α₁M+α₄C+K(u_i+1)；

S6：Calculate remainder：

R=(α₁(u_i+1-u_i)-α₂v_i-α₃a_i)M+C(α₄〔u_i+1-u_i〕+α₅v_i+α₆a_i)+R(u_i+1)-f_i+1；

Judge r square whether less than TOL, if it is, the current u of output_i+1、v_i+1And a_i+1；Otherwise, execution step S7；

S7：Solve equation A (xu_i+1)=r, then according to formula below updates u_i+1、v_i+1And a_i+1：

Wherein,

Perform step S2.

Preferably, step D specifically includes step：

D1：According to the displacement on 8 summits in the hexahedral element, and initial coordinate p before deformation '₀, p ’₁... p '₇, obtain 8 current coordinates in summit：p₀, p₁... p₇；Wherein, p '₀Coordinate be (x, y, z)；

D2：Calculate a=| P '-p '₀| x/Q, b=| P '-p '₀| y/Q, c=| P '-p '₀|z/Q；Wherein, Q represents described hexahedro The length of side of body unit；P ' represents the initial coordinate on the summit of the leaf model intermediate cam grid before deformation；

D3：Interpolation weights f on 8 summits of the hexahedral element are calculated respectively₀, f₁... f₇：

f₀=(1-a) × (1-b) × (1-c)；

f₁=a × (1-b) × (1-c)；

f₂=a × b × (1-c)；

f₃=(1-a) × b × (1-c)；

f₄=(1-a) × (1-b) × c；

f₅=a × (1-b) × c；

f₆=a × b × (1-c)；

f₇=(1-a) × b × c；

D4：The current coordinate P in the summit of the leaf model intermediate cam grid is calculated according to the following equation：

(3) beneficial effect

The plant leaf blade wilting analogy method based on physical model of the present invention, by threedimensional model voxel collection, base are constructed The Mass-spring Model built in the voxel collection is multilamellar, is more pressed close to the internal structure of real plants blade, and matter Connect more stable between point spring；Employing implicit expression Newmark integral algorithm can be simulated under big step-length, it is ensured that Model can carry out real-time visual Simulation；Simultaneously by can be with real time modifying damping parameter, stiffness factor with user mutual Deng attribute, the Mass-spring Model for meeting blade physical arrangement feature can be more accurately constructed, satisfaction carries out plant leaf blade and moves The needs of state deformation simulative.

Description of the drawings

Fig. 1 is the flow chart of the plant leaf blade wilting analogy method described in the embodiment of the present invention based on physical model；

Fig. 2 is the voxel collection schematic diagram of leaf model described in the embodiment of the present invention；

Fig. 3 is the spring structure schematic diagram included in single hexahedron；

Fig. 4 is the structural representation of the Mass-spring Model；

Fig. 5 a are the first design sketchs when the cucumber leaves that simulation is obtained are wilted；

Fig. 5 b are the second design sketchs when the cucumber leaves that simulation is obtained are wilted.

Specific embodiment

With reference to the accompanying drawings and examples, the specific embodiment of the present invention is described in further detail.Hereinafter implement Example is not limited to the scope of the present invention for illustrating the present invention.

Fig. 1 be described in the embodiment of the present invention based on physical model plant leaf blade wilting analogy method flow chart, such as Fig. 1 Shown, methods described includes step：

A：The voxel collection of three-dimensional plant leaf model is built, the voxel collection includes multiple for simulating plant blade mesophyll The hexahedral element of internal cell.

Step A specifically includes step：

A1：The segmentation of Octree space is carried out to the triangle gridding in leaf model and obtains space Octree.

A2：The spatial surrounding box of the leaf model is calculated, the spatial surrounding box is divided into into multiple hexahedral elements.

A3：All hexahedral elements are traveled through, judges each described hexahedral element with the space Octree successively Subtree space whether intersect, if it is, retaining the current hexahedral element；Otherwise, the current hexahedron is given up Unit；

A4：The traversal institute hexahedral element with a grain of salt, judges successively each described hexahedral element and the subtree Whether the triangle gridding in space intersects, if it is, choosing the current hexahedral element；Otherwise, give up current described Hexahedral element.

A5：The voxel collection of the leaf model is built using all hexahedral elements chosen, and records all The index information on the summit of the triangle gridding included in the hexahedral element chosen, so as to by the summit of triangle gridding and phase The hexahedral element answered is corresponding, is that subsequent treatment carries out data preparation.Fig. 2 is the body of leaf model described in the embodiment of the present invention Element collection schematic diagram, as shown in Fig. 2 the voxel is concentrated includes multiple hexahedral elements, and the hexahedral element is Square.

B：The voxel collection is converted into into Mass-spring Model.The Mass-spring Model is made up of multiple springs.Fig. 3 is The spring structure schematic diagram included in single hexahedron, as shown in figure 3, the spring by the hexahedral element in any two Individual summit connects and composes, and the summit of the both ends of the spring is denoted as particle.Fig. 4 is the structural representation of the Mass-spring Model, As shown in figure 4, the Mass-spring Model is on the basis of the voxel collection being made up of the hexahedral element, by the leaf model The spring of multiple interconnections is divide into, by the design of this spring structure, matter in the leaf model can be effectively improved The stability of point spring connection.

The internal force formula of the spring is as follows：

E=k_s(|x_ij|-L)²；

Wherein, E represents the energy between particle i and particle j, x_ijThe vector difference of particle i and particle j is represented, L represents particle Initial distance between i and particle j, k_sRepresent the stiffness factor of the spring comprising particle i and particle j.

The physical kinetics equation of the Mass-spring Model is as follows：

Ma+D (u, v)+R (u)=F；

Wherein, u ∈ R³ⁿRepresent the displacement of particle；v∈R³ⁿFor the first derivative of u, the speed of particle is represented；a∈R³ⁿFor u Second dervative, represent particle acceleration；M∈R^3n×3nRepresent the mass matrix of particle；D (u, v) ∈ R³ⁿRepresent that particle is subject to Damping force；R(u)∈R³ⁿRepresent the internal force that particle is subject to；R represents real number space；N represents particle in the Mass-spring Model Number；F∈R³ⁿRepresent the external force that particle is subject to.

The computing formula of internal force R (u) that the particle is subject to is as follows：

Wherein, x_iRepresent the three dimensional space coordinate of particle i.

The computing formula of damping force D (u, v) that the particle is subject to is as follows：

D (u, v)=(α M+ ω K (u)) v；

Wherein, α and ω represent the non-negative damping parameter of the damping characteristic for controlling the particle spring system, can basis Practical situation real time modifying.Preferably, α be 0, ω be 1.K(u)∈R^3n×3nTangent stability matrix is represented, computing formula is as follows：

Wherein, x_jThe three dimensional space coordinate of particle j is represented, I represents unit matrix.

C：The Mass-spring Model is applied after external force, simulation calculates the position on each summit in the hexahedral element Move.

Step C specifically includes step：

C1：Initial setting up is carried out to the stiffness factor of each spring in the Mass-spring Model.For belonging to vein position Spring, it is 0.5～1.0 that stiffness factor arranges scope, from petiole more close to device of spring stiffness coefficient it is bigger, and the spring at edge Stiffness factor is less.For the spring of mesophyll position, it is 0.3～0.5 that stiffness factor arranges scope, random distribution.For each Hexahedral element, can be subject to a downward gravity, and when wilting, the summit of hexahedral element is subject to one to unit center Contractility (power represents blade during dehydration wilting, turgescence change caused by the contraction of vacuole), while hexahedral element Comprising spring stiffness factor diminish (excursion of the stiffness factor of the spring of vein position be 0.1～0.2, mesophyll The excursion of the stiffness factor of the spring of position be 0.2～0.4).With the reduction of water content, the amplitude that stiffness factor diminishes Become big, the summit of hexahedral element contractility and become big.

C2：The Mass-spring Model is applied after external force, using implicit expression Newmark integration method to the physical kinetics Equation carries out numerical integration, is calculated the position on each summit (particle in i.e. described spring model) in the hexahedral element Move.By the way that preferably convergence effect can be reached using the implicit expression Newmark integration method, it is allowed to which user adopts the larger time Step-length is simulated.

The implicit expression Newmark integration method is as follows：

Input includes：Displacement components u of the particle at the i moment_i, speed vi and acceleration a_i, the i+1 moment be subject to it is outer Power f_i+1, iterationses j_max, and iteration error upper limit TOL, time step t；

Output includes：In the displacement components u at i+1 moment_i+1, speed v_i+1And acceleration a_i+1；

And the implicit expression Newmark integration method includes step：

S1：Setting u_i+1Initial value be u_i, the initial value of variable j is 0；

S2：Whether j judges j less than j from after Jia 1_max, if it is, execution step S3, otherwise, exports current u_i+1、 v_i+1And a_i+1；

S3：According to current u_i+1, calculate the internal force R (u that particle is subject to_i+1), and tangent stability matrix K (u_i+1)；

S4：Calculate Local Damping Matrix C=(α M+ ω K (u))；

S5：Calculate symmetrical matrix A=α₁M+α₄C+K(u_i+1)；

S6：Calculate remainder：

R=(α₁(u_i+1-u_i)-α₂v_i-α₃a_i)M+C(α₄〔u_i+1-u_i〕+α₅v_i+α₆a_i)+R(u_i+1)-f_i+1；

Judge r square whether less than TOL, if it is, the current u of output_i+1、v_i+1And a_i+1；Otherwise, execution step S7；

S7：Solve equation A (xu_i+1)=r, then according to formula below updates u_i+1、v_i+1And a_i+1：

Wherein,

Perform step S2.

The displacement of each particle in the spring model is obtained by the implicit expression Newmark integration method, institute has also just been obtained State the displacement on each summit in hexahedral element.

D：According to the displacement on each summit in the hexahedral element, interpolation goes out the top of the leaf model intermediate cam grid The coordinate of point.Each described hexahedral element includes 8 summits.

Step D specifically includes step：

D1：According to the displacement on 8 summits in the hexahedral element, and initial coordinate p before deformation '₀, p ’₁... p '₇, obtain 8 current coordinates in summit：p₀, p₁... p₇；Wherein, p '₀Coordinate be (x, y, z).

D2：Calculate a=| P '-p '₀| x/Q, b=| P '-p '₀| y/Q, c=| P '-p '₀|z/Q；Wherein, Q represents described hexahedro The length of side of body unit；P ' represents the initial coordinate on the summit of the leaf model intermediate cam grid before deformation.

D3：Interpolation weights f on 8 summits of the hexahedral element are calculated respectively₀, f₁... f₇：

f₀=(1-a) × (1-b) × (1-c)；

f₁=a × (1-b) × (1-c)；

f₂=a × b × (1-c)；

f₃=(1-a) × b × (1-c)；

f₄=(1-a) × (1-b) × c；

f₅=a × (1-b) × c；

f₆=a × b × (1-c)；

f₇=(1-a) × b × c.

D4：The current coordinate P in the summit of the leaf model intermediate cam grid is calculated according to the following equation：

According to the changing coordinates on the summit of each triangle gridding of leaf model, plant leaf blade wilting just can be simulated Form afterwards.Fig. 5 a are the first design sketchs when the cucumber leaves that simulation is obtained are wilted；Fig. 5 b are the cucumber leaves that simulation is obtained The second design sketch during wilting.As shown in figure 5 a and 5b, by the method for the invention, can be with more life-like simulation plant Form when blade is wilted.

Plant leaf blade wilting analogy method described in the embodiment of the present invention based on physical model, by threedimensional model tectosome Element collection, the Mass-spring Model built based on the voxel collection is multilamellar, is more pressed close to the internal structure of real plants blade, And connect between mass spring more stable；Employing implicit expression Newmark integral algorithm can be simulated under big step-length, Ensure that model can carry out real-time visual Simulation；Simultaneously by can be with real time modifying damping parameter, strength with user mutual The attributes such as degree coefficient, can more accurately construct the Mass-spring Model for meeting blade physical arrangement feature, and satisfaction carries out plant The needs of blade dynamic deformation simulation.

Embodiment of above is merely to illustrate the present invention, and not limitation of the present invention, about the common of technical field Technical staff, without departing from the spirit and scope of the present invention, can also make a variety of changes and modification, therefore all The technical scheme of equivalent falls within scope of the invention, and the scope of patent protection of the present invention should be defined by the claims.

Claims (9)

1. a kind of plant leaf blade wilting analogy method based on physical model, it is characterised in that including step：

A：Build three-dimensional plant leaf model voxel collection, the voxel collection include it is multiple for simulating plant blade mesophyll inside The hexahedral element of cell；

B：The voxel collection is converted into into Mass-spring Model；

C：The Mass-spring Model is applied after external force, simulation calculates the displacement on each summit in the hexahedral element；

D：According to the displacement on each summit in the hexahedral element, interpolation goes out the summit of the leaf model intermediate cam grid Coordinate；

Step A specifically includes step：

A1：The segmentation of Octree space is carried out to the triangle gridding in leaf model and obtains space Octree；

A2：The spatial surrounding box of the leaf model is calculated, the spatial surrounding box is divided into into multiple hexahedral elements；

A3：All hexahedral elements are traveled through, the son of each described hexahedral element and the space Octree is judged successively Whether tree space intersects, if it is, retaining the current hexahedral element；Otherwise, the current hexahedral element is given up；

A4：The traversal institute hexahedral element with a grain of salt, judges successively each described hexahedral element with the subtree space Whether interior triangle gridding intersects, if it is, choosing the current hexahedral element；Otherwise, give up current described hexahedro Body unit；

A5：The voxel collection of the leaf model is built using all hexahedral elements chosen, and records all choosing The hexahedral element in the index information on the summit of triangle gridding that includes.

2. the method for claim 1, it is characterised in that the Mass-spring Model is made up of multiple springs；The bullet Spring is connected and composed by any two summit in the hexahedral element, and the summit of the both ends of the spring is denoted as particle.

3. method as claimed in claim 2, it is characterised in that the internal force formula of the spring is as follows：

E=k_s(|x_ij|-L)²；

Wherein, E represents the energy between particle i and particle j, x_ijRepresent particle i and particle j vector difference, L represent particle i and Initial distance between particle j, k_sRepresent the stiffness factor of the spring comprising particle i and particle j.

4. method as claimed in claim 3, it is characterised in that the physical kinetics equation of the Mass-spring Model is as follows：

Ma+D (u, v)+R (u)=F；

Wherein, u ∈ R³ⁿRepresent the displacement of particle；v∈R³ⁿFor the first derivative of u, the speed of particle is represented；a∈R³ⁿFor the two of u Order derivative, represents the acceleration of particle；M∈R^3n×3nRepresent the mass matrix of particle；D(u,v)∈R³ⁿWhat expression particle was subject to hinders Buddhist nun's power；R(u)∈R³ⁿRepresent the internal force that particle is subject to；R represents real number space；N represents particle in the Mass-spring Model Number；F∈R³ⁿRepresent the external force that particle is subject to.

5. method as claimed in claim 4, it is characterised in that the computing formula of internal force R (u) that the particle is subject to is as follows：

$R\left(u\right)=\frac{\∂E}{\∂x_i}=k_s\left(\right|x_{ij}|-L)\frac{x_{ij}}{\left|x_{ij}\right|}=k_{s}u\frac{x_{ij}}{\left|x_{ij}\right|};$

Wherein, x_iRepresent the three dimensional space coordinate of particle i.

6. method as claimed in claim 4, it is characterised in that the computing formula of damping force D (u, v) that the particle is subject to is such as Under：

D (u, v)=(α M+ ω K (u)) v；

Wherein, α and ω represent the damping parameter of the damping characteristic for controlling the particle spring system；Represent just Stability matrix is cut, computing formula is as follows：

$K\left(u\right)=\frac{\∂R\left(u\right)}{\∂x_j}=k_s\frac{x_{ij}{x_{ij}}^T}{{x_{ij}}^T{x}_{ij}}+k_s(1-\frac{L}{\left|x_{ij}\right|})(I-\frac{x_{ij}{x_{ij}}^T}{{x_{ij}}^T{x}_{ij}});$

Wherein, x_jThe three dimensional space coordinate of particle j is represented, I represents unit matrix.

7. method as claimed in claim 6, it is characterised in that step C specifically includes step：

C1：Initial setting up is carried out to the stiffness factor of each spring in the Mass-spring Model；

C2：The Mass-spring Model is applied after external force, using implicit expression Newmark integration method to the physical kinetics equation Numerical integration is carried out, the displacement on each summit in the hexahedral element is calculated.

8. method as claimed in claim 7, it is characterised in that the implicit expression Newmark integration method is specific as follows：

Input includes：Displacement components u of the particle at the i moment_i, speed v_iAnd acceleration a_i, in the external force that the i+1 moment is subject to f_i+1, iterationses j_max, and iteration error upper limit TOL, time step t；

Output includes：In the displacement components u at i+1 moment_i+1, speed v_i+1And acceleration a_i+1；

And the implicit expression Newmark integration method includes step：

S1：Setting u_i+1Initial value be u_i, the initial value of variable j is 0；

S2：Whether j judges j less than j from after Jia 1_max, if it is, execution step S3, otherwise, exports current u_i+1、v_i+1With And a_i+1；

S3：According to current u_i+1, calculate the internal force R (u that particle is subject to_i+1), and tangent stability matrix K (u_i+1)；

S4：Calculate Local Damping Matrix C=(α M+ ω K (u))；

S5：Calculate symmetrical matrix A=α₁M+α₄C+K(u_i+1)；

S6：Calculate remainder：

R=(α₁(u_i+1-u_i)-α₂v_i-α₃a_i)M+C(α₄(u_i+1-u_i)+α₅v_i+α₆a_i)+R(u_i+1)-f_i+1；

Judge r square whether less than TOL, if it is, the current u of output_i+1、v_i+1And a_i+1；Otherwise, execution step S7；

S7：Solve equation A (xu_i+1)=r, then according to formula below updates u_i+1、v_i+1And a_i+1：

$\left\{\begin{array}{c}{u}_{i+1}=u_i+{\mathrm{xu}}_{i+1}\\ v_{i+1}={\mathrm{\α}}_4(u_{i+1}-u_i)+{\mathrm{\α}}_5{v}_i+{\mathrm{\α}}_6{a}_i\\ a_{i+1}={\mathrm{\α}}_1(u_{i+1}-u_i)-{\mathrm{\α}}_2{v}_i+{\mathrm{\α}}_3{a}_i\end{array}\right.;$

Wherein,

Perform step S2.

9. the method for claim 1, it is characterised in that step D specifically includes step：

D1：According to the displacement on 8 summits in the hexahedral element, and initial coordinate p before deformation '₀,p’₁,……p’₇, Obtain 8 current coordinates in summit：p₀,p₁,……p₇；Wherein, p '₀Coordinate be (x, y, z)；

D2：Calculate a=| P '-p '₀| x/Q, b=| P '-p '₀| y/Q, c=| P '-p '₀|z/Q；Wherein, Q represents the hexahedron list The length of side of unit；P ' represents the initial coordinate on the summit of the leaf model intermediate cam grid before deformation；

D3：Interpolation weights f on 8 summits of the hexahedral element are calculated respectively₀,f₁,……f₇：

f₀=(1-a) × (1-b) × (1-c)；

f₁=a × (1-b) × (1-c)；

f₂=a × b × (1-c)；

f₃=(1-a) × b × (1-c)；

f₄=(1-a) × (1-b) × c；

f₅=a × (1-b) × c；

f₆=a × b × (1-c)；

f₇=(1-a) × b × c；

D4：The current coordinate P in the summit of the leaf model intermediate cam grid is calculated according to the following equation：

$P=\underset{i=0}{\overset{7}{\Σ}}{p}_i{f}_i.$