An Isogeometric Cloth Simulation Based on Fast Projection Method

Xuan Pengand Chao Zheng

1Department of Mechanical Engineering,Suzhou University of Science and Technology,Suzhou,215009,China

2Simshow High Technology Co.,Ltd.,Chengdu,610042,China

3Department of Mechanical and Industrial Engineering,The University of Iowa,Iowa City,IA 52242-1527,USA

ABSTRACT A novel continuum-based fast projection scheme is proposed for cloth simulation.Cloth geometry is described by NURBS,and the dynamic response is modeled by a displacement-only Kirchhoff-Love shell element formulated directly on NURBS geometry.The fast projection method,which solves strain limiting as a constrained Lagrange problem,is extended to the continuum version.Numerical examples are studied to demonstrate the performance of the current scheme.The proposed approach can be applied to grids of arbitrary topology and can eliminate unrealistic over-stretching efficiently if compared to spring-based methodologies.

KEYWORDS Cloth simulation;isogeometric analysis;strain limiting;fast projection

1 Introduction

Cloth simulation has many practical applications,such as computer-aided garment design,character animation,and electronic e-commerce.Terzopoulos et al.[1]were among the first to develop a physical model for use in the simulation of cloth.Breen et al.[2–4]developed a spring-based model,and their initial motivation was to accurately represent complex fabric behaviors using nonlinear springs.Provot [5] modeled cloth by employing linear springs characterized by low stiffness and he achieved acceptable results with high efficiency.Eischen et al.[6] modeled cloth using finiteelement shell theory.Baraff et al.[7] proposed to employ implicit time integration to increase the time step,while the Newton iteration was suggested to be halted at the first step to achieve high efficiency.Choi et al.[8] proposed an immediate buckling model to ensure that the Jacobi matrix in the implicit method was definite.Bridson et al.employed a velocity Verlet algorithm[9]to increase the time increment.In addition to the cloth model,the collision response was also largely improved.Baraff et al.[7]simulated collision as a velocity constraint and Bridson et al.[10]proposed a decoupled bullet-proof collision scheme.

The garment industry is now beginning to use virtual simulation for prototyping[11,12].However,as reported in[13],there are still many challenges in the application of virtual simulation.One of the challenges is to develop a high-efficiency continuum approach.The continuum approach has many advantages over the spring-based approach.For example,the material properties of the continuum approach are independent of the topology of the grid.Another attraction of the continuum approach is that many studies of material and geometric nonlinearity have been conducted in this field.However,conventional finite-element shell theory has low efficiency due to two reasons:(1)it needs more degrees of freedom for the same grid,and(2)the collision response is more complex.

Another challenge of cloth simulation is how to efficiently enforce realistic strain on cloth.One of the characteristics of fabric is that the bending stiffness is far lower than the in-plane stiffness.A consequence of this property is that in-plane deformation of practical cloth in most cases is negligible if compared to the out-of-plane deformation.However,using high physical in-plane stiffness introduces significant difficulty in simulation.For explicit methods,higher in-plane stiffness requires smaller time increments.In Barraf et al.[7] semi-implicit method,the frictional damping is proportional to the in-plane stiffness.Thus,to reduce the computational burden,in many practical cloth simulations,the value of the in-plane stiffness is artificially reduced to enable the simulation of more complex problems in reasonable time.However,this approach can introduce unrealistic overstretching of the cloth.

To overcome this undesired side effect,Provot[5]proposed an explicit method of strain limiting that restores the over-stretched springs by adjusting the particle position directly.Because adjusting the position of one spring may result in over-stretching of another spring,an iteration is required.Both Jacobi and Gauss-Seidel iterations[10]are utilized,but neither one can guarantee convergence.Implicit methods of strain limiting have also been proposed [4,14].These approaches are mainly based on the constrained Lagrange method and consider in-plane strain as a constraint condition.Goldenthal et al.[15] proposed a fast projection method that can solve the constrained Lagrange problem much more efficiently.The implicit method requires solving a linear system for each iteration,but it can converge very quickly.The fast projection method in [15] constrains the spring length of a weft and warp aligned quadrilateral grid,but this kind of grid is difficult to obtain for practical garments.To obtain a grid-independent strain-limiting scheme,one must resort to a continuum model.Studies of continuum-based strain limiting are limited at present.Thomaszewski et al.[16]proposed explicit strain limiting for triangular elements.

Isogeometric analysis has been proposed to bridge computer-aided design (CAD) and analysis seamlessly[17].The IGA owns the salient features such as higher-order continuity and exact geometry preservation etc.,which provides an effective solution for problems that conventional finite element methods are not proper qualified to solve[18–20],thus it has received much attention in many fields.For cloth simulation,a systematic method is proposed by the authors [21].This method uses the rotational-free Kirchhoff-Love shell model [22,23],wherein CAD geometry is directly utilized in analysis.Recent developments in cloth-like simulation involve large deformation shells or membrane modeling [24–26].The present work proposes a continuum version of a strain-limiting scheme by applying the fast projection method,such that the shell model for cloth simulation can be solved quickly with an acceptable resultant shape.Another advantage of this model is that it performs the simulation directly on the NURBS surface,which is widely implemented in CAD software.This property ensures the high fidelity of the simulated cloth and provides convenience of interactive design in three-dimensional(3D)space.

The rest of this paper is organized as follows.Section 2 briefly reviews the rotation-free Kirchhoff-Love Shell on a NURBS basis.Section 3 outlines the continuum-based fast projection method for(trimmed) NURBS geometry.Finally,four example problems are solved to validate the proposed approach and their results are presented in Section 4.Conclusions are presented in Section 5.

2 NURBS Kirchhoff-Love Shell

Kirchhoff-Love shell theory assumes the following:

· The normal to the undeformed middle surface remains straight and perpendicular to the deformed middle surface.

· The transverse normal stress is small compared with other normal stress components and may be neglected.

· The thickness of the shell is small compared to the other dimensions.

· The displacements of any given point on the shell are small in comparison to the thickness.

These assumptions are a good approximation for fabrics in which the energetic contribution from transverse shear is negligibly small compared to the bending and in-plane energy.Hence,the kinetics are completely characterized by the surface strain and curvature,which are determined by the surface geometry.For numerical computation,the theory can lead to a displacement-only formulation,which does not involve rotational degrees of freedom,thus increasing the efficiency of the method.A Kirchhoff-Love shell element is not commonly used in traditional finite-element analysis because constructing aC1continuous surface is difficult for some element topologies.However,with NURBS geometry,it is straightforward to construct surfaces ofC1or even higher-order continuity.For this and other considerations,the NURBS Kirchkoff-Love shell is utilized for cloth modeling.

In the present NURBS Kirchkoff-Love shell element,the primary unknowns are the displacements of the control points.No rotational degree-of-freedoms are introduced.

2.1 Kinematics

The NURBS formulation below follows that of Kiendl et al.[22].We use the same set of NURBS basis functions to parameterize the reference and current configurations of the cloth surface:

Here,the QIs are the reference coordinates of control points and the qIare the current coordinates.TheNIs are the basis functions of the NURBS surface.The knot parametersξ=(ξ1,ξ2)serve as the convected coordinates whereby a fixed pairξrepresents the same material point throughout deformation.These two coordinates induce two convected surface basis vectors a1=x,ξ1and a2=x,ξ2spanning the tangent plane at every point of the surface.A line element in the current configuration is thus represented asdx=a1dξ1+a2dξ2.In the reference configuration,the surface bases are denoted{A1,A2},and a line element is given bydX=A1dξ1+A2dξ2.The basis vectors are illustrated in Fig.1.

The unit normal a3(in the current configuration)is

With respect to the convected basis vectors,the surface deformation tensor C and Green-Lagrangian strain E take the forms

respectively.The surface curvature tensorκis defined by

It is convenient to use a local ortho-normal basis to perform the element computations presented later.To this end,we introduce a pair of orthonormal basespoint-wise in the tangent plane spanned by{A1,A2}.With respect to the new bases,we express a line element asdX=Locally,(d¯X1,d¯X2)are related to(dξ1,dξ2)via

Derivatives of a basis functionNin local physical coordinates follow the chain rule:

In the current configuration,the basesconvected from the physical basisare

With respect to the physical basis,the Green-Lagrangian strain assumes the form-δαβ).The physical components of the curvature tensor can be obtained by the transformation

2.2 Constitutive Law

Cloth response is typically inelastic,exhibiting anisotropic properties and a small to moderate amount of hysteresis [27–29].Since the focus of the present work is on geometry,we simplify the constitutive description by using an isotropic elastic model.The in-plane strain of fabrics is usually small (＜2%);however,since large rotation is involved,the use of finite strain is necessary.Because of the small strain range,any mechanically sound finite-strain constitutive model should reasonably describe the in-plane response.Here,we use a linear anisotropic relation [7] between the (in-plane)Piola-Kirchhoff stress S and Green-Lagrangian strain E:

in whichIcan be replaced by“weft,”“warp,”or“shear,”andkIis a material constant.

For the bending model,we employ a nonlinear bending function[29]for curvatures in weft and warp directions,

whereIcan be replaced by“weft”or“warp,”andB0andκ0are material constants.The bending curve is depicted in Fig.2.The shape mimics the ascending portion of a typical fabric bending in a Kawabata test[30].The bending moment,in general,depends on curvature change;here,the reference curvature is taken to be zero.

Figure 2:Moment-curvature curve in principal space

2.3 Element Equation

External forces acting on a piece of cloth normally include a body forceρb and damping force fd(per unit surface area);that is,in general,a function of the cloth velocity,and traction forces ¯tprescribed on the boundary edge.Following the textile community convention,we use the surface densityρto describe the mass distribution;ρ=ρ0h,whereρ0is the 3D density andhis cloth thickness.The weak form of the dynamics equilibrium equation is given by

In the NURBS representation,and thus

where I is the 3×3 identity matrix.The variation of Green-Lagrangian strain,in Voigt form,δE=(δE11,δE22,2δE12),is derived as

where

Similarly,from the definition of curvature,we can derive

where

In the above,j=||a1×a2||is the area stretch,and(a1,a2)are the bases dual to(a1,a2)satisfying aα·aβ=δαβand aα·a3=0.

Substituting Eqs.(11),(12),and(14)into Eq.(10)yields the discrete dynamic equation

where

For low-speed air drag,we assume that fd=-η.,whereηis a viscosity constant;in this case,

2.4 Time Integration

Our strain-limiting scheme is independent of time integration.Here,we use the velocity Verlet scheme that was first applied to cloth simulation by [9].The flow process for time integration is as follows:

1.Predict average velocity and candidate configuration attn+1,

3 Continuum-Based Fast Projection Method

3.1 Constrained Lagrange Method

The fast projection method begins with the constrained Lagrange problem.For the general coordinates r=(qT,λT),the constraint Lagrange potential is given by

in which C is a set of constraint conditions that stops in-plane stretching.

The Euler-Lagrange equation is:Thus,we have,

Supposing that we check the constraint condition at the end of a time step,C(qn+1=0),and use the velocity Verlet scheme for internal dynamics,we then have:

Now,we use the splitter:

1.Predict candidate configurationby internal dynamics(without constraints)

This sub-step can be replaced by any time-integration scheme.

2.Correct the candidate configuration by qn+1=+Δq,and solveΔq from

Eqs.(25a)and(25b)are equivalent toδW=0,where

The second term ofWprojects q into the constraint manifold,while the first term minimizes the change of kinematic energy in this projecting process.All together,q will be projected to the“closest”point on the constraint manifold.

Linearizing Eqs.(25a)and(25b)by the Newton method,we have the following flow process:

Δq[0]=0,λ[0]=0,for j=0,1,2,...,do

1.SolveΔq anddλfrom

2.Correct q andλby

The iteration exits when the C is smaller than a given tolerance.We note here that we must solve a(3n+m)-dimensional asymmetric linear system for each iteration.

3.2 Fast Projection Scheme

We note that Eq.(25a) keeps the property of “closest” and Eq.(25b) keeps the constraint conditions.In the actual problem,we have a high requirement on constraint conditions but a low requirement on the property of“closest.”Goldenthal et al.[15]suggested solving Eq.(25a)using the Euler method while solving Eq.(25b)using the Newton method.We expand C by a Taylor series at q[j]and ignore the quadratic and higher-order terms

Substituting Eq.(29)into Eq.(25a)yields

Eq.(25b)is still linearized by the Newton method

Substituting Eq.(30)into Eq.(31)yields the basic equation of the fast projection method

The iterative fast projection process flow is as follows:

Δq[0]=0,for j=0,1,2,...

1.Solvedλfrom

2.Correct q by

The iteration exits when the C is smaller than the given tolerance.The advantage of the fast projection method is that:(1)we only need to solve an m-dimensional S.P.D.linear system for each iteration,and(2)it converges faster.

3.3 Spring-Based Fast Projection

For the spring-mass method,the length of each spring is a constraint,C=[C1,C2,...],and

in whichliandLiare the current and referent lengths of theithspring,respectively.For triangular mesh,because the total number of springs is approximately 3 times the total number of nodes,the deformation of the entire mesh is locked and even curvature is not allowed.Thus,quad-dominated mesh is required.

3.4 Continuum-Based Fast Projection

For the continuum-based approach,the first task is to select sampling points.We tried to select the Gaussian points as sampling points,but that does not work well.Because the number of cells are close to the number of control points,and if there are three or more constraint conditions on each cell,the model will be locked.Thus,2×2 or more Gaussian points are not allowed.We then tried one Gauss point and the average strain,but both schemes display the hourglass problem.By adding a certain term to the averaged strain scheme,we might have been able to eliminate the hourglass problem,but instead we decided to try a different approach.We selected vertices of a knot mesh as sampling points,as shown in Fig.3.For trimmed NURBS,we constrain the intersection points between the trimming curve and knot mesh as well.The benefit of this scheme is that it does not have locking or hourglass problems,and it can minimize the bandwidth of the coefficient matrix

The constraint condition C is expressed as CT=[CT1,CT2...],in which Ciis a constraint on one sampling point.For each sampling point,we can set,at most,two constraints.Considering that the shear stiffness of fabric is far lower than the stretch stiffness in weft and warp directions,we constrain the strain in weft and wart directions,given by:

in whichEweftandEwarpare the strains in the weft and warp directions of the fabric,respectively.Noting that the local basesandare defined as any pair of orthogonal vectors in the tangent plane spanned by A1and A2,a simple way to computeEweftandEwarpis to setandaligned with the weft and warp directions.

Figure 3:Constraint points

3.5 Stress Reconstruction

The constraint energy termλ·C,in essence,is the nodal integration form of strain energy.Thus,the stressshould include two terms:the stress from strain S(E)and that from constraint force.The stress on grid pointiafter fast projection is given by

whereaiis the nodal area.

However,we found a checkerboard pattern in the values ofλ.This is mainly due to the lack of domain integration.Introducing the integration version of constraint energy may resolve the problem,but we believe that is to costly as the bandwidth will be heavily enlarged.Therefore,we introduced a simple way to reconstruct the stress: by an interpolation,which makes it easy to obtain the stress on Gaussian points.Then,by an extrapolation from Gaussian points,we can obtain the nodal stress without a checkerboard problem.

4 Examples

4.1 Corner Kidnapped Cloth

A piece of cloth in the x-y plane is subject to constraints at two corners and will swing under gravity in the z direction.The cloth is represented by a second-order NURBS patch with 100 control points.The bending parameters areB0=3.3×10-3N·m,κ0=30 m-1.The mass density isρ=0.117 kg/m2and the damping constant isη=0.351 kg/(m2·s).Fabric thickness ish=0.001177 m.For the in-plane model,we assumek=kweft=kwarp=2.5kshearand test three cases:(a)k/ρ=100,strain limit=0.01;(b)k/ρ=5000 without strain limiting;and(c)k/ρ=100 without strain limiting.

The simulation results of different Young’s moduli att=0.6 s(the cloth passes the vertical plane for the first time near this time) are shown in Fig.4.The simulation time and strain summary are shown in Table 1.It is observed that case (c) has obvious unrealistic stretching.The results of cases(a)and(b)are close,but the time increment of case(a)is much smaller than that of case(b).Because there is no contact involved,the time savings of a large time increment is not very obvious;however,strain limiting still obtains a speedup of 2.5 times.

Table 1:Summary of corner kidnapped cloth simulation

Figure 4:Corner kidnapped cloth at t=0.6 for cases(a)–(c)

4.2 Draping of Soft Armor

The draping of a soft armor was simulated in this example.The armor is represented by a secondorder NURBS patch with 1,616 control points.The initial configuration is obtained by virtual try-on simulation and is shown in Fig.5.The bending parameters areB0=1×10-4N·m,κ0=30 m-1.The mass density isρ=0.1177 kg/m2and the damping constant isη=0.353 kg/(m2·s).Fabric thickness ish=0.001177 m.The friction coefficient is assumed to beμ=0.2.For the in-plane model,we assumek=kweft=kwarp=2.5kshearand create three cases:(a)k/ρ=50,strain limit=0.01;(b)k/ρ=5000 without strain limiting;and(c)k/ρ=50 without strain limiting.

Figure 5:Initial configuration of soft armor

To more clearly show the results,the simulation results of upper-and lower-body armor att=0.5sare shown in Figs.6 and 7,respectively.The simulation time and strain summary are shown in Table 2.For the upper-body armor,case (c) obtains an unacceptable shape,while for lower-body armor all three schemes obtain an acceptable shape.By checking the maximum and average strain,it is found that both cases (a) and (b) obtain a low level of strain,but the time cost of case (b) is much greater than case (a).This demonstrates that the strain-limiting method can reduce strain level by reducing extra CPU cost rather than by pure high stiffness.

Table 2:Summary of soft armor draping simulation

Figure 6:Upper-body armor after draping simulation for cases(a)–(c)

Figure 7:Lower-body armor after draping simulation for cases(a)–(c)

4.3 Draping of a Skirt

This example simulates the draping process of a skirt.The initial configuration of the skirt is obtained by a try-on simulation and is shown in Fig.8.The entire model contains 960 control points.The woman’s body is represented by a discrete mesh of 17,068 cells.The bending parameters areB0=1.0×10-5N·m andκ0=30 m-1.The mass density isρ=0.118 kg/m2and the damping constant isη=0.354 kg/(m2·s).The fabric thickness ish=0.001177 m.A friction coefficientμ=0.8 is assumed.For the in-plane model,we assumek=kweft=kwarp=2.5kshearand create three cases: (a)k/ρ=50,strain limit=0.01;(b)k/ρ=3000 without strain limiting;and(c)k/ρ=50 without strain limiting.

The garment shape att=0.8 s,when the vibrations are mostly damped out,is shown in Fig.9.The simulation time and strain summary are shown in Table 3.It is observed that low stiffness without strain limiting obtains an unacceptable shape.WhenE/ρ=3000 and strain limiting is turned off,the simulation shape looks better,but local over-stretching still occurs.We did not try a higher stiffness because it takes too much time.However,by using the proposed strain-limiting scheme,both the simulation shape and strain values satisfied requirements while keeping the time cost very low.

Table 3:Summary of skirt draping simulation

4.4 Patch Test

Finally,we conducted a standard patch test to check the stress obtained from the fast projection method.A 1 m-by-1 m square was fixed on the top edge and we applied 1 N/m of uniform force on the bottom edge.As the problem itself is based on dynamic analysis,we modeled the static patch test by applying a large damping factor and waiting until the vibration was damped out.Fig.10 shows the static stress obtained from analysis with or without the fast projection method.It is observed that for both cases the stresses are close to 1 Pa,which means that the fast projection method cannot only eliminate unrealistic strain,but can also provide a reliable stress field.

Figure 8:Initial configuration of skirt

Figure 9:Skirt after draping simulation for cases(a)–(c)

Figure 10:Patch test stress

5 Conclusions

In this work,a rotational-free Kirchhoff-Love shell-based isogeometric analysis was outlined for cloth simulation.To overcome the numerical burden caused by high in-plane stiffness,a continuum version of the fast projection method was applied.The highlights of this work are the following:

· Compared with spring-based models,the constraint directions of strain limiting are independent of grid lines.This implies that the fast projection method can be applied to a grid with arbitrary topology.

· Examples show that the present scheme can eliminate unrealistic over-stretching efficiently,while the stress field remains reliable.

· The trimmed NURBS patches are used directly for geometry,indicating seamless application of CAD and analysis.

Future work will focus on developments of efficient numerical methods to handle localized features,e.g.,wrinkles,to facilitate a“real-time”simulation of cloth.

Funding Statement:Chao Zheng thanks the support from Sichuan Science and Technology Program[Grant No.2021JDRC0007].

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.