Int J Mech Mater Des (2005) 2:225–243 DOI 10.1007/s10999-006-9004-0

A new skeletal model for fabric drapes K. Y. Sze Æ X. H. Liu

Received: 5 January 2006 / Accepted: 16 May 2006 / Published online: 13 July 2006
Springer Science+Business Media B.V. 2006

Abstract Fabric drapes are typical large displacement, large rotation and small strain problems. Compared with continuum shell finite element methods, methods that skeletonize the fabric sheet into a set of interconnected nodes appear to be more popular in drape problems with extensive wrinkles. These skeletal methods may resort to particle mechanics and formulate the elastic energy or the equations of motion in terms of the node-to-node distance and the angles between the straight lines joining adjacent nodes. Alternatively, beam elements can be employed to skeletonize the fabric sheet at the expense of using rotational in addition to translational nodal DOFs. In this paper, a new skeletal model based on the small strain and curvature assumptions is devised. In contrast to most, if not all, skeletal models for fabric drape simulation, all the stretching, shearing and bending energies in the model are simple polynomials of the gridpoint displacement. Fabric drape problems with

K. Y. Sze (&) Æ X. H. Liu Department of Mechanical Engineering, The University of Hong Kong, Pokfulam Road, Hong Kong SAR, P.R. China e-mail: [email protected] Present Address: X. H. Liu Department of Mechanics, Huazhong University of Science & Technology, Wuhan 430074, P.R. China

extensive wrinkles are examined. The presence of sharp fold, seam and cut in the undeformed fabric sheet as well as fabric-to-solid contact are considered. The predicted appearances are pleasant and conform to real life observations. Keywords Beam Æ Geometric nonlinearity Æ Grid Æ Fabric drape

Introduction There is an ever-growing desire to automate the design and manufacture of clothing and the related products. Images produced by drape simulation software can be used for demonstrating the garment to designers. The reduction in time and effort at the iterative design stage can be substantial. Furthermore, the notion of a customer selecting a garment for purchase via internet, showing the garment draped over his/ her body model, has been suggested as a potential e-commerce marketing method. To this end, a key element in related technologies is an efficient and robust computational method which can effectively predict the mechanics response of fabrics. In the last two decades, considerable research efforts have been made in fabric drape simulation. Majority of the work can broadly be categorized

123

226

Int J Mech Mater Des (2005) 2:225–243

‘‘- - -’’ connecting a node and its six second closest neighborhoods, see Fig. 2b. Figure 3 depicts yet another way of connecting the nodes. The gridlines are again aligned with either the warp or the weft. Associated with the typical node ‘‘O’’, the stretching energy is stored by four linear springs along the gridlines and the shearing energy is stored by four diagonal linear springs, see Fig. 3b, as in Fig. 1d. The bending energy is stored by two linear springs connected alternating nodes along the gridlines and another two linear springs connecting alternating nodes diagonal to the grid, see Fig. 3c (Choi and Ko 2002). If the bending energy associated with the twisting curvature defined with respect to the gridlines is neglected, the two diagonal springs can be excluded (Lin et al. 2003; Villard and Borouchaki 2005). Same as Zhang and Yuen (2001), Choi and Ko (2002), Lin et al. (2003), and Villard and Borouchaki (2005), only linear springs are used in the work of Meyer et al. (2001). On the other hand, Stylios et al. (1996) developed a node-bar model based on a physical analogue to a deep shell system and Ascough et al. (1996) skeletonized a fabric sheet into a network of thin-beam finite elements. Teng et al. (1999) proposed a model which provides a systematic way for reducing fabric sheets into orthogonal grids. Some notions of the finite volume method were invoked in the derivation of the model. Same as the conventional interacting particle model, torsional springs are used to store the shearing and bending energies. By adopting a corotating frame and assuming that both the strain and curvature are small, the authors have recently been able to express the stretching and bending energies as simpler functions of the nodal displacements (Sze 2002; Sze and Liu 2003). Unlike skeletal models, fabric surfaces appear explicitly in the continuum-based methods.

into the geometrical (Hinds et al. 1991; Rodel et al. 1998; Dai et al. 2001) and mechanistic methods (Breen et al. 1992; Gan et al. 1995; Kang and Yu 1995; Chen and Govindaraj 1995; House et al. 1996; Eberhardt et al. 1996; Stylios et al. 1996; Ascough et al. 1996; Eischen et al. 1996; Tan et al. 1999; Teng et al. 1999; Zhang and Yuen 2001; Meyer et al. 2001; Gong et al. 2001; Choi and Ko 2002; Sze 2002; Sze and Liu 2003; Lin et al. 2003; Villard and Borouchaki 2005; Etzmuss et al. 2003). The geometrical methods are mainly adopted by the computer graphics community and focus on animating a cloth-like appearance whereas mechanistic methods resort to mechanistic principles such as particle mechanics, energy methods and shell theories. Furthermore, mechanistic methods can be grouped into the skeletal and the continuum methods. In the skeletal methods, a fabric sheet is skeletonized into a network of interconnected nodes or grid-points where kinetic d.o.f.s are defined. A well-known example is the conventional interacting particle model (Breen et al. 1992; House et al. 1996; Eberhardt et al. 1996) in which the gridlines are orthogonal and are aligned with either the warp or the weft, see Fig. 1. Associated with the typical particle ‘‘O’’, the stretching, shearing and bending energies are stored by the four linear, four torsional and two torsional springs as shown in Fig. 1b, c and e, respectively. It is possible to replace the four torsional springs for storing the shear energy with four diagonal linear springs, see Fig. 1d. Figure 2 portrays another way of aligning the nodes and each of which is connected to 12 neighborhoods (Zhang and Yuen 2001). The stretching and shear energies are stored by the linear springs ‘‘—’’ connecting a node and its six closest neighborhoods whereas the bending energy is stored by the linear springs

(a)

(b)

(c)

(d)

(e)

Fig. 1 (a) A typical particle ‘‘O’’ connected to its four neighborhoods, (b) linear springs for storing the stretching energy, (c) torsional and (d) diagonal linear springs for storing the shearing, and (e) torsional springs for storing the bending energy

123

Int J Mech Mater Des (2005) 2:225–243

227

Fig. 2 (a) The grid pattern. (b) A typical particle ‘‘O’’ connected to its twelve neighborhoods: — are linear springs storing stretching and shearing energy, - - - are linear springs storing the bending energy

The normal Green strain which is a quadratic function of the nodal displacement is used to quantifying the stretching energy. By adopting some un-invoked small strain and small curvature approximations, the bending and shear energies can be expressed in terms of simply polynomial functions of the nodal displacements. Computational trials indicate that the prediction of the previous and the present models are very close in all the examples considered in Sze and Liu (2003), yet the latter is not only simpler but also more efficient than the former. In this paper, new examples with extensive wrinkles are presented with sharp fold, seam and contact constraints. The predicted appearances appear to be realistic.

Methods of solution for nonlinear elastic systems

(a)

(b)

(c)

Fig. 3 (a) A typical particle ‘‘O’’ connected to its 16 neighborhoods, (b) linear springs for storing the stretching and shearing energies, (c) linear springs for storing the bending energy

Among the continuum-based methods, most of them employed shell finite element models (Gan et al. 1995; Kang and Yu 1995; Chen and Govindaraj 1995; Eischen et al. 1996; Tan et al. 1999; Gong et al. 2001; Etzmuss et al. 2003). These include thin triangular plate/shell elements in which the transverse shear energy is ignored (Tan et al. 1999; Gong et al. 2001; Etzmuss et al. 2003) and the bending energy may simply be computed in terms of the angle of non-coplanarity between adjacent elements. Others are degenerated shell elements (Gan et al. 1995; Kang and Yu 1995; Chen and Govindaraj 1995; Eischen et al. 1996) among which advanced finite element techniques including assumed natural strain method and reduced integration with stabilization control can be noted. To the best knowledge of the authors, previous attempts by shell finite element models are rather restricted to drape problems with limited wrinkles. In this paper, a new skeletal model will be devised. Unlike the previously proposed model (Sze 2002; Sze and Liu 2003), the present model does not involve transformation between the corotational and the global coordinate frames.

Consider an elastic deformable body whose elastic energy E can be expressed as a function of the displacements of the nodes distributed over the body, i.e. E ¼ EðUÞ

ð1Þ

where U ¼ fUi g is the system nodal displacement vector. In general, the total potential energy of the system can be expressed as: P ¼ EðUÞ
PT U

ð2Þ

in which P ¼ fPi g is the system external force vector that accounts for the prescribed external forces. The static equilibrium condition can be obtained by minimizing the potential energy, i.e. @P ¼ @U




@P @Ui

 ¼0

or; equivalently; F
P ¼ 0 ð3Þ

where F ¼ fFi g ¼ f@E=@Ui g

is system internal force

vector arising from the internal forces. To solve the above equation, the full Newton– Raphson iterative scheme can be used, i.e.

123

228

Int J Mech Mater Des (2005) 2:225–243

Uð0Þ ¼ 0; and

DUðiÞ ¼ ðKjU¼UðiÞ Þ
1 ðP
FjU¼UðiÞ Þ

Uðiþ1Þ ¼ UðiÞ þ DUðiÞ

ð4Þ

in which UðiÞ is the ith iterative solution of U; DUðiÞ is the iterative refinement of UðiÞ ; K ¼ ½Kij  ¼ ½@Fi =@Uj  ¼ @FT =@U ¼ @=@Uð@E=@UÞT is the system tangential stiffness matrix. The iteration can be stopped when the iterative displacement refinement DU or the unbalanced force F
P is sufficiently small. Alternatively, one may also use a minimization algorithm such as the conjugate gradient method to obtain the equilibrium solution. The major advantage of these algorithms is that only the total potential energy P and its gradient @P=@U ¼ F
P have to be known whereas the Hessian or system tangential stiffness matrix is not required.

Integration zone and elastic energies under consideration Similar to the many skeletal models (Breen et al. 1992; House et al. 1996; Eberhardt et al. 1996; Teng et al. 1999; Dai et al. 2001; Meyer et al. 2001; Choi and Ko 2002; Sze 2002; Sze and Liu 2003; Lin et al. 2003; Villard and Borouchaki 2005), the gridlines in the present model are orthogonal and are aligned with the warp and weft directions. The translational dofs are defined at the nodes or grid-points. With the use of the finite volume concept (Teng et al. 1999), node-based integration zones are defined and a typical one for node ‘‘O’’ away from the boundary is shown in Fig. 4a (Sze and Liu 2003). The elastic energy to be considered in the node-based integration zone can be expressed as: O O O O EO ¼ EO AB
e þ EAB
j þ ECD
e þ ECD
j þ EBD
c O O þ EO DA
c þ EAC
c þ ECB
c

123

ð5Þ

Fig. 4 (a) Node-based integration zone for a typical node O. Node-based integration zones for (b) stretching and bending energies along warp, (c) stretching and bending energies along weft and (d) shearing energy. h denotes node-based integration zone, O denotes node and – – – denote gridline O where EO AB
e and EAB
j are stretching and bending energies arising from, respectively, the tensile strain and the curvature along A-O-B; O EO CD
e and ECD
j are stretching and bending energies arising from, respectively, the tensile O strain and curvature along C-O-D; EO BD
c , EDA
c , O O EAC
c and ECB
c are the shearing energies in the first, second, third and fourth quadrants of the integration zone, respectively. It can be noted that the transverse shear energy, the coupling energy effect of the tensile strains along the warp and weft directions, the coupling energy effect of the curvatures along the warp and weft directions as well as the energy associated with the twisting curvature are neglected for simplicity and computational efficiency (Breen et al. 1992; House et al. 1996; Eberhardt et al. 1996; Teng et al. 1999; Dai et al. 2001; Sze and Liu 2003; Lin et al. 2003; Villard and Borouchaki 2005). For the nodes away from the boundary, the edges defining an integration zone are parallel and equidistant from their adjacent gridlines, see Fig. 4b–d. For the nodes on and next to the boundary, the integration zones for the stretching/bending energy along the warp direction, the stretching/bending energy along the weft

Int J Mech Mater Des (2005) 2:225–243

direction and the shearing energy are different as depicted in Fig. 4b–d, respectively.

229

w ¼ NA wA þ NO wO þ NB wB ¼ NO wO

ð6Þ

where Bending energy As reviewed in Sect. 1, there are different ways of quantifying the bending energy in the skeletal models. With reference to fabric thread AOB which is deformed into A¢O¢B¢ in Fig. 5, one may use the difference between \AOB and \A0 O0 B0 (see Fig. 1 and Breen et al. 1992; House et al. 1996; Eberhardt et al. 1996; Teng et al. 1999; Dai et al. 2001) or the difference between the lengths of AB and A¢B¢ (see Figs. 2, 3, and Zhang and Yuen 2001; Choi and Ko 2002; Lin et al. 2003; Villard and Borouchaki 2005). Though Zhang and Yuen (2001), Choi and Ko (2002), Lin et al. (2003), and Villard and Borouchaki (2005) do not mention why the change in length is preferred rather than the change of angle, the authors’ experience indicate that the latter involves anti-trigonometric function. The resulting bending energy and its derivative(s) required by the method of solution are computationally less efficient. In the previous work of the authors, a corotational frame is defined by A¢, O¢ and B¢ with A¢ being the origin as shown in Fig. 5. The first coordinate of the frame passes through A¢, Oc and B¢ which are colinear. In particular, O¢ and Oc are coincident if A¢O¢B¢ is strain free. The interpolated transverse deflection w with respect to the corotational frame is:

Fig. 5 AOB is deformed to A¢O¢B¢ and s is the corotating coordinate which is also the physical length measured from A along AOB or from A¢ along A¢OcB¢

ðs
aÞðs
lÞ ; al sðs
aÞ NB ¼ : bl

NA ¼

NO ¼


ðs
lÞs ; ab

Using the small curvature assumption, the radius of curvature would be: 1 d2 w 2 ¼
2 ¼ wO : R ds ab

ð7Þ

Thus, there are at least two ‘‘non-angular’’ methods to quantify the bending energy. With reference to Fig. 5, they are uB and wO. To further examine the two methods, the fabric thread in Fig. 6 is considered. AOB is bent into a circular arc A¢O¢B¢ of radius R. It can be seen that u and w are employed in the two methods for quantifying the bending energy. By assuming that the stretching strain of thread is small, one can write R ¼ l; u þ 2R sin

 ¼ l; 2

w þ R cos

 ¼ R: 2 ð8Þ

Under the small curvature assumption, l/R > 1. By Taylor’s expansion, the expressions in (8) lead to:   l l 1 l3 l
u ¼ 2R sin ’ 2R
and 2R 2R 6 8R3   l 1 l2 ’R 1
R
w ¼ R cos 2R 2 4R2

Fig. 6 An extensible fabric strip draped into a circular arc

123

230

Int J Mech Mater Des (2005) 2:225–243

or 1 24 ’ 3u R2 l

and

1 8 ¼ w R l2

ð9Þ 2

For positive u, both u and w can be used to quantify the bending energy which is proportional to 1/R2. However, when the fabric is pulled, u is negative and the bending energy being proportional to 1/R2 would also be negative which is physically impossible and can induce computational problems. Most, if not all, of the works based on the first method employ u2 to quantify the bending energy. As u2 ~ 1/R4, such a choice is also not sensible from mechanistic point of view and should be restricted to animation but not simulation. For any straight segment of a fabric thread, the transverse displacement component does not enter the bending energy. As a result, additional device should be required to subdue the relevant instability. Though (7) is a very simple and concise expression, wO can only be obtained by going through the transformation required in the corotational formulation. Further simplifications will be induced here to facilitate the determination of the bending energy. In fabric drape, the external load is the gravity which, in general, will not induce much stretching strain. With reference to local displacement components u and w in Fig. 5, one should have w ? u. Thus, wO ’ jOc O0 j. As the position vectors of O and Oc can be expressed as: b a b a XO ¼ XA þ XB and XOc ¼ XA0 þ XB0 l l l l b a ¼ ðUA þ XA Þ þ ðUB þ XB Þ; l l one gets b a w ’ jOc O0 j ¼ jXO0
XOc j ¼ jUO
UA
UB j l l ¼ jNUAOB
j j: ð10Þ

123

In the last expressions, N ¼ ½
bl I3 ; I3 ;
al I3  and UAOB
j ¼ fUTA ; UTO ; UTB gT . The bending energy is: DjAB L EO ¼ AB
j

aþb=2 Z

2 2

j2 ds¼

12LDjAB 2 a2 b2

a=2 aþb=2 Z

w2 ds¼

12lLDjAB T UAOB
j NT NUAOB
j 2 a2 b2

a=2

ð11Þ where DjAB is the flexural rigidity along A-B or the warp. Differentiation yields @EO 2lLDjAB T AB
j ¼ N NUAOB
j and @UAOB
j a2 b2  T @EO 2lLDjAB T @ AB
j ¼ N N @UAOB
j @UAOB
j a2 b2

ð12Þ

which will be assembled to the form the system internal force vector and tangential stiffness matrix. The internal force vector and the tangential stiffness matrix arising from the bending energy EO CD
j along CD or the weft can be derived in a similar manner. Stretching energy In our previous work, the stretching strain along A¢O¢B¢ is measured by using the displacement with respect to the corotating frame attached to A¢, O¢ and B¢ as shown in Fig. 5. Using the displacement measured with respect to the corotation frame, the Green strain was employed, i.e. e¼

    du 1 du 2 1 dw 2 þ þ ds 2 ds 2 ds

ð13Þ

where s is a physical length measured from A. The third local displacement component normal to the plane formed by A¢, O¢ and B¢ is zero by definition. In this paper, the general Green strain is used to exempt the coordinate transformation required in the corotational formulation, i.e.

Int J Mech Mater Des (2005) 2:225–243

1 2¼ X;Ts U;s þ U;Ts U;s : 2

231

ð14Þ

The position vector X and the displacement U = {U, V, WT} along AOB can be interpolated as: X ¼ NA XA þ NO XO þ NB XB ¼ fs; 0; 0gT ; U ¼ fNi Ui ; Ni Vi ; Ni Wi g

T

ð15Þ

where Ns have been defined in (6) and repeated subscripts imply summation over A, O and B. For instance, Ni Ui ¼ NA UA þ NO UO þ NB UB . Substitution of (15) into (14) gives, 1 e ¼ N 0i Ui0 þ ðN 0i N 0j Ui Uj þ N 0i N 0j Vi Vj þ N 0i N 0j Wi Wj Þ 2 1 ð16Þ ¼ N 0i Ui þ N 0i N 0j ðUi Uj þ Vi Vj þ Wi Wj Þ 2 in which N¢i denotes the derivative of Ni with respect to s. The stretching energy would be DeAB L EO AB
e ¼

aþb=2 Z

2 2

e2 ds ¼

1 LDjAB 2 2

a=2



aþb=2 Z

½N 0i N 0j Ui Uj

a=2

þN 0i N 0j N 0h ðUi Uj þ Vi Vj þ Wi Wj ÞUh 1 þ N 0i N 0j N 0h N 0k ðUi Uj þ Vi Vj þ Wi Wj Þ 4  ðUh Uk þ Vh Vk þ Wh Wk Þds: where DeAB is the tensile rigidity along A-B or the warp. The above energy can be re-written as 1 1 EO AB
e ¼ Mij Ui Uj þ Mijh ðUi Uj þ Vi Vj þ Wi Wj ÞUh 2 2 1 þ Mijhk ðUi Uj þ Vi Vj þ Wi Wj ÞðUh Uk 8 þ Vh Vk þ Wh Wk Þ ð17Þ

It is trivial that the above terms are symmetry with respect any pair of indices and there are 6, 10 and 15 independent Mijs, Mijhs and Mijhks, respectively. By invoking the symmetry, differentiation yields @EO AB
e ¼Mmi Ui @Um 1 þ Mmij ð3Ui Uj þVi Vj þWi Wj Þ 2 1 þ Mmijh ðUi Uj þVi Vj þWi Wj ÞUh ; 2 @EO AB
e ¼Mmij Ui Vj @Vm 1 þ Mmijh ðUi Uj þVi Vj þWi Wj ÞVh ; 2 @EO AB
e ¼Mmij Wi Uj @Wm 1 þ Mmijh ðUi Uj þVi Vj þWi Wj ÞWh ; 2 @ 2 EO AB
e ¼Mmn þ3Mmni Ui @Um @Un 1 þ Mmnij ð3Ui Uj þVi Vj þWi Wj Þ; 2 @ 2 EO AB
e ¼Mmni Ui @Vm @Vn 1 þ Mmnij ðUi Uj þ3Vi Vj þWi Wj Þ; 2 @ 2 EO AB
e ¼Mmni Ui @Wm @Wn 1 þ Mmnij ðUi Uj þVi Vj þ3Wi Wj Þ; 2 @ 2 EO AB
e ¼Mmni Vi þMmnij Ui Vj ; @Um @Vn @ 2 EO AB
e ¼Mmni Wi þMmnij Wi Uj ; @Wm @Un @ 2 EO AB
e ¼Mmijh Vi Wj : ð18Þ @Vm @Wn

in which ðMij ; Mijh ; Mijhk Þ ¼

LDjAB 2

aþb=2 Z

a=2

ðN 0i N 0j ; N 0i N 0j N 0h ; N 0i N 0j N 0h N 0k Þ ds:

By grouping the above derivatives, the following internal force vector and stiffness matrix arising from the stretching energy EO AB
e can be formed:

123

232

Int J Mech Mater Des (2005) 2:225–243

9 8 u v w u  1 > = < Mu þ 2 ðM u þ M v þ M wÞ þ M u þ M u > ; ¼ Mv u þ M v > @UAOB
e > ; : Mw u þ M w 2 3 M þ 3Mu þ Muu þ M Mv þ Muv Mw þ Mwu  T O @ @EAB
e 6 7 ¼4 Mu þ Mvv þ M Mvw 5 @UAOB
e @UAOB
e u ww  symmetry M þM þM @EO AB
e

ð19Þ

in which

u ¼ fUA ; UO ; UB gT ; v ¼ fVA ; VO ; VB gT ; w ¼ fWA ; WO ; WB gT ; UAOB
e ¼ fuT ; vT ; wT gT ; ðMÞij ¼ Mij ;

ðMu Þij ¼ Mijh Uh ;

ðMvv Þij ¼ Mijhk Vh Vk ;

ðMv Þij ¼ Mijh Vh ;

ðMww Þij ¼ Mijhk Wh Wk ;

ðMw Þij ¼ Mijh Wh ;

ðMuv Þij ¼ Mijhk Uh Vk ;

ðMuu Þij ¼ Mijhk Uh Uk ; ðMvw Þij ¼ Mijhk Vh Wk ;

1 ðMwu Þij ¼ Mijhk Wh Uk and M ¼ ðMuu þ Mvv þ Mww Þ: 2

The expressions in (19) can be assembled to form the system internal force vector and system tangential stiffness matrix. The internal force vector and the tangential stiffness matrix arising from the stretching energy EO CD
e along CD or the weft can be derived in a similar manner.

Examples on two-dimensional drapes In applying the minimum potential energy or virtual work principle to continua, displacement continuity is a pre-requisite. Though the presented procedure avoids the use of rotational d.o.f.s which require special treatment in large displacement/rotation analysis, displacement continuity is not secured as the interpolation zone is larger than the integration zone. Hence, the procedure should only be considered as an improvement of the particle-based methods. But still, it is important that the present procedure can yield predictions consistent with other more rigorous approaches such as finite element method. The derivation in the previous section has been sufficient for testing the present model in one-dimensional fabric and beam bending

123

problems. To this end, only a single row of nodes will be considered and all nodes will be quipped with the translational d.o.f.s along the X- and Zdirections but not the Y-direction. Numerical solutions presented in this section are obtained by the conventional full Newton– Raphson iteration method, see (4). The following force convergence criterion is employed : jF
Pj=jPj610
4 :

ð20Þ

To mimic the clamped boundary condition, the slave node technique commonly used in the finite difference method is employed. The slave and the master nodes are symmetrically placed about the clamped node, see Fig. 7a. Displacement of the former node is constrained so that its position and the position of the master node are symmetric with respect to the un-deformed transverse normal at the clamped node. A clamped fabric strip deforms under its own weight This problem has been considered by Kang and Yu (1995), Teng et al. (1999), Sze and Liu (2003).

Int J Mech Mater Des (2005) 2:225–243

233

0 x+ ++ 0+ x x +

1

2

x (cm)

+x+x

-1

x x+ +

x

+x +x

(a)

+

-2

z (cm)

+x

+x

4

5

weft

(b)

Kang & Yu (exptl.) Teng et al previous, 12 nodes previous, 22 nodes present, 12 nodes present, 22 nodes

+x

+x

-3

3

xx+

warp +x

+x

-4 +x

+x

-5

(c)

Fig. 7 (a) Strip modeled by 12 nodes, ‘‘•’’ denotes the master ‘‘m’’ and slave ‘‘s’’ nodes. (b) strip modeled by 22 nodes and (c) the deformed strips under gravity

The fabric properties are listed in Table 1. Strips of unit width cut along the warp and weft directions are clamped at one end and the freehanging length is 5 cm. The 5 cm length is modeled by 11 and 21 non-uniformly distributed nodes as depicted in Fig. 7a, b. An extra slave node is added for implementing the clamped boundary condition. Gravitational pull is applied in a single load step and the deformed strips are depicted in Fig. 7c. The predictions of the previous and present models are close. They are also in good agreement with the experimental results. Predictions yielded by more nodes in the previous and present models are graphically indistinguishable from the ones obtained by using 22 nodes. The results of Teng et al. (1999) were obtained by using 56 non-uniformly distributed nodes. Cantilever subjected to end shear A 10 m long, 1 m wide and 0.1 m thick cantilever is subjected to a concentrated end shear force with its maximum magnitude equal to

Table 1 Mechanical properties of the fabric under consideration (Kang and Yu 1995)

15 kN, see Fig. 8a. The cantilever is made of an isotropic material whose elastic modulus and Poisson’s ratio are 1.2 GPa and 0, respectively. To assess the reliability of the present results, the problem is analyzed by the commercial finite element code ABAQUS using 6 · 60 elements (ABAQUS 2000; Sze and Liu 2003). The fournode CPS4I incompatible plane element is chosen due to its excellent bending response. In the present method, the 10 m length is modeled by 11 and 21 uniformly distributed nodes. Again, one slave node is required for implementing the clamped condition. With the full load applied in 30 equal load increments, Fig. 8b plots the end shear force against the free-end displacements. More elements and nodes have been employed but the yielded predictions are graphically indistinguishable from the ones obtained by 6 · 60 elements and 22 nodes, respectively. The converged predictions, with respect to the number of nodes or elements solutions, of ABAQUS, the previous and the present methods are very close. When 12 nodes are used, the prediction of present method is closer to the

Tensile rigidity, D
(gf/cm)

Bending rigidity, Dj (gf cm2/cm)

Warp

Weft

Warp

Weft

1118.2

759.5

0.083

0.063

Shear rigidity, Dc (gf/cm)

Weight (gf/cm2)

Thickness, h(cm)

41.8

0.019

0.0593

123

234

Int J Mech Mater Des (2005) 2:225–243 16

x

14

+

12

ABAQUS, 6x60 elements previous, 12 nodes

x +

x+

x +

previous, 22 nodes present, 12 nodes present, 22 nodes

x+

x+

x+

x+

x+

Load P (kN)

x+

x+

x+

10

x+

x+

x+

x+

8

x+

U x+

W x+

x+

6

x+

x+ x+

4

(a)

x+ x+

x+

x+

x+

2

+ x

0 + x 0

+ x

x+ x+

+ x 1

2

3

4

5

6

7

8

9

End Deflections (m)

(b) Fig. 8 (a) A cantilever loaded by end shear force and (b) the load–deflection curves

converged solutions than that of the previous method. Shearing energy The next energy component to be considered is shearing energy. To this end, the shearing energy in the first quadrant of the integration zone in Fig. 1b is considered and it can be written as: 1 bdDc 2 cBD 2 4

EO BD
c ¼

ð21Þ

c

where D is the shear rigidity and cBD is the shear strain. With the small strain assumption, Breen et al. (1992), House et al. (1996), Eberhardt et al. (1996), Teng et al. (1999), and Dai et al. (2001) calculate the shear strain using the shear angle, i.e. p XTB0 O0 XTD0 O0 cBD ¼
cos
1 ; 2 jXB0 O0 j  jXD0 O0 j Dc bd 2 c EO BD
c ¼ 2 4 BD

ð22Þ

where XPQ = XP – XQ, X¢P = XP+UP is the displaced position of node P and UP is the displacement of node P. For higher computational efficiency and in compliance with the small strain assumption, the following approximation was adopted in our previous work (Sze and Liu 2003):

123

cBD ’ sin cBD ¼

XTB0 O0 XTD0 O0 : jXB0 O0 j  jXD0 O0 j

ð23Þ

In this paper and in line with the small strain assumption, the above expression is further approximated as: cBD ’

XTB0 O0 XTD0 O0 XT 0 0 XT 0 0 XT 0 0 XT 0 0 ’ BO DO ¼ BO DO : bd jXB0 O0 j  jXD0 O0 j jXBO j  jXDO j ð24Þ

From (21) and (24), it can be derived that @EO BD
c @UBOD

¼

bdDc @cBD cBD ; 4 @UBOD

T bdDc  @c  @c T @  @EO BD
c BD BD ¼ @UBOD @UBOD 4 @UBOD @UBOD @  @cBD T  þcBD ; @UBOD @UBOD ð25Þ where 9 9 8 8 X D0 O0 > > > = @c = < UB > < 1 BD ¼
XB0 O0
XD0 O0 ; UBOD ¼ UO ; > > > ; @UBOD bd > ; : : UD XB0 O0 2 3 0
I3 I3  T @ @cBD 16 7 ¼ 4
I3 2I3
I3 5: @UBOD @UBOD bd I3
I3 0

Int J Mech Mater Des (2005) 2:225–243

235

The expressions in (25) can be assembled to form the system internal force vector and system tangential stiffness matrix. The internal force vector and the tangential stiffness matrix arising O from the shearing energies EO DA
c , EAC
c and O ECB
c for other quadrants of the integration zone can be derived in a similar manner. With all the energies derived, the present model can be applied to drape simulation of fabric sheets. Three-dimensional drapes In this section, square fabric draped over square table, square fabrics draped over circular tables and square fabrics draped over spheres are considered. The fabric properties listed in Table 1 are adopted. Fabric sheet draped over a square table In this example, a 20 · 20 cm fabric sheet is draped over a 10 · 10 cm square table. Owing to symmetry, only a quarter of the fabric sheet is modeled. The simulation is attempted by the previous and the present skeletal models with 11 · 11, 21 · 21 and 51 · 51 uniformly distributed nodes. The nodes in contact with the table are fully restrained and no slave node is used. Kang and Yu (1995) also attempted the same problem with the four-node MITC4 degenerated shell element and 17 · 17 uniformly distributed

nodes was employed to cover a quarter of the fabric sheet. The predicted top views of the draped sheet together with the experimental result are shown in Fig. 9a. The predictions yielded by the previous and the present skeletal models are indistinguishable regardless of the nodal density. It can be seen that the predictions yielded by 51 · 51 nodes match quite well with the experimental result which is, however, slightly nonsymmetric. More nodes have also been employed but the predictions are indistinguishable from that by 51 · 51 nodes. Considerable difference can be noted between the finite element prediction and experimental result. The former, however, may not have converged. Figure 9b shows the top view of the fabric sheet predicted by 31 · 31 nodes with smaller nodal spacing at the vicinity of the table. The prediction is obtained by the present model but it is graphically indistinguishable from that of the previous model (Sze and Liu 2003; Teng et al. 1999). Fabric sheet draped over a circular table As a matter of fact, all the 16 three-dimensional drape examples which include square sheets over circular tables, circular sheets over square tables and circular sheets over circular tables considered by the previous model in Sze and Liu (2003), have been attempted by the present one. All these examples involve extensive wrinkling. Despite of the fact that several

(a) 10

(b)

Y (cm)

5

previous & present, 11x11 nodes previous & present, 21x21 nodes previous & present, 51x51 nodes

0

Kang & Yu (exptl.) Kang & Yu (FEM 17x17 nodes)

-5

-10 -10

-5

0

5

10

X (cm)

Fig. 9 Top views of a square fabric sheet draped over a square table obtained by (a) uniformly nodal distribution and (b) 31 · 31 uniformly distributed nodes

123

236

Int J Mech Mater Des (2005) 2:225–243

fully restrained. To better present the restraining effect of the circular table, the simple scheme of adding a smaller number of fully restrained nodes to the intersecting points of the table boundary and the gridlines devised is adopted, see (Sze and Liu 2003) for details. Figure 10 shows the top and isometric views predicted by the previous and present models. The readers can probably see that the predictions by the two models are exceptional close and the difference can only be noted after careful examination.

additional simplifications based on small strain and small curvature assumptions have been incorporated in the present model, the predictions of the two models are surprisingly close. In this light, they are not repeated here. To illustrate how close the model predictions can be, a new example which involves more extensive wrinkling than all the reported ones is considered. Here, a 320 · 320 cm fabric sheet is draped over a Ø160 cm table. The sheet is essentially modeled by 81 · 81 uniformly distributed nodes. The nodes over the tables are 100

100

Y

Z

Y

X

Z

50

X

50

Y (cm)

Y (cm)

0

-50

0

-50

-100 -100

-50

0

50

100

-100 -100

-50

(a) predicted by the previous model

50

100

(b) predicted by present model Z

Z Y

0

X (cm)

X (cm)

X

(c) predicted by the previous model

Y

X

(d) predicted by present model

Fig. 10 The predicted top and isometric views for a 320 · 320 cm square fabric sheet draped over a Ø160 cm table by the previous and the present models

123

Int J Mech Mater Des (2005) 2:225–243

237

Fabric sheet draped over a sphere

Ci ðUÞ¼

Figure 11 shows a 180 · 180 cm wool fabric sheet draped symmetrically over a Ø80 cm sphere. 91 · 91 nodes are employed in modeling the fabric sheet. With the centres of the sheet and the sphere taken to be respectively the coordinate origin and (0,0,R), the contact constraints for problem are: X

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðXi þUi Þ2 þðYi þVi Þ2 þðZi þWi
RÞ2
R ð26Þ T

for i = 1,...,8281 where {Xi, Yi, Zi} is the undeformed nodal position vector and {Ui, Vi, Wi}T is the nodal displacement vector. The constraints can be enforced by penalty method and the total potential energy in (2) is augmenting to be:

Y

P ¼ EðUÞ
PT U þ a

Z

N X

½½Ci 2

ð27Þ

i¼1

in which a denotes a large positive penalty factor whereas ½½ ¼  for \0 and ½½ ¼ 0 for  > 0. The penalized minimization problem can be solved by either Newton–Raphson’s method or a minimization algorithm such as the conjugate gradient method. Figure 12 shows the top, side

Fig. 11 A 180 · 180 cm2 wool fabric sheet draped over a Ø80 cm sphere. 91 · 91 nodes are employed

Y

X

X

Z

Y Y

X Z Z

Fig. 12 Top, side and isometric views of a square wool fabric draped over a Ø80 cm sphere

123

238

Int J Mech Mater Des (2005) 2:225–243

The reinforcements are 6 cm wide and form a cross. The tensile rigidity, shear rigidity, bending rigidity, and weight of the reinforcements are 5, 5, 20, and 2 times of those of fabric, respectively. The computed top, side and isometric views are shown in Fig. 14. Compared with other parts of the fabric sheets, the stiffness of the two reinforcements uphold themselves from falling down.

and isometric views of the computed solutions. It can be seen that the draped shape is symmetric about the X- and Y-axes. Fabric sheet with reinforced boundary draped over a sphere The fabric sheet considered in Sect. 8.3 is reinforced along its four edges. The reinforcement is 3 cm wide. Its tensile rigidity, shear rigidity, bending rigidity and weight are respectively 5, 5, 20 and 2 times that of the plain fabric sheet. The computed top, side and isometric views are shown in Fig. 13. Comparing with Fig. 12, the symmetry of the draped shape about the X- and Y-axes has been lost due to the presence of the reinforcement.

Fabric sheet with a cut In this example, the fabric sheet considered in Sect. 8.3 is cut along the warp direction. The cut is 8 cm from the boundary and 90 cm in length, see Fig. 15a. The computed top, side and isometric views are shown in Fig. 15b–d.

Fabric sheet with reinforced axes of symmetry draped over a sphere

Two fabric sheets jointed by a seam

In this example, the fabric sheet considered in Sect. 8.3 is reinforced along its axes of symmetry.

In this example, two 92 · 180 cm fabric sheets are jointed by a 2 cm seam as shown in Fig. 16. Y

X

Z

X

Y

Y

X

Z Z

Fig. 13 Top, side and isometric views of fabric sheet with 3 cm reinforced boundary

123

Int J Mech Mater Des (2005) 2:225–243

239 Y

X

Z

X X

Y

Y Z Z

Fig. 14 Top, side and isometric views of a fabric sheet with 6 cm wide reinforcements along its axes of symmetry

Each of two sheets is modeled by 47 · 91 nodes and no integration zone crosses the boundary of the two sheets. Over the seam, the nodes are stitched and each pair of the stitched nodes share the same nodal dofs. In the resting state, the jointed fabric sheets are hanging vertically downwards side by side, see Fig. 16a. All the non-stitched nodes are then raised to the same vertical level as the lower end of the seam and the nodal spacing is kept at 2 cm, see Fig. 16b. Before the drape starts, the whole fabric sheet is vertically shifted until it is marginally in contact with the Ø80 cm sphere. Figure 16c–f depict the top, side and isometric views of the final draped shape. Fabric Sheet with a fold In this example, a fold along the weft direction in the middle of the 180 · 180 cm fabric sheet is considered. Each side of the fold is modeled as

a 90 · 180 cm sheet with no integration zone crossing the fold and 46 · 91 nodes are employed. To simulate the fold, the top rows of nodes in the two sheets are collapsed and the next two rows of nodes are jointed by linear springs with zero unstretched lengths, see Fig. 17a. The energy in the linear spring joining node I and node J is: k EIJ ¼ ðXI þ UI
XJ
UJ ÞT 2  ðXI þ UI
XJ
UJ Þ

ð28Þ

in which k is the spring constant and is taken to be 2 gf/cm. In the rest state, the two sides of the fold are hanging downwards side by side. Thus, XI = XJ and  k k I
I3 EIJ ¼ ðUI
UJ ÞT ðUI
UJ Þ ¼ UTIJ 3 UIJ
I3 I3 2 2 ð29Þ

123

240

Int J Mech Mater Des (2005) 2:225–243

(a)

Y

(b) X

(c)

Z

(d)

X

Y

X

Z

Y

Z

Fig. 15 A fabric sheet with a 90 cm cut

where UIJ ¼ ½UTI ; UTJ T . Differentiations give:  I3
I3 dEIJ ¼k UIJ ; dUIJ
I3 I3    I3
I3 d dEIJ T ¼k dUIJ dUIJ
I3 I3

ð30Þ

which can be assembled to form the system internal force vector and tangential stiffness matrix. With the nodal spacing maintained at 2 cm, all nodes are raised to the same vertical level and the centre of the fold is marginally in contact with the upper pole of the Ø80 cm sphere, see Fig. 17b. In other words, the fold is unfolded and only the linear springs, which represent the restoring stiffness of the fold, are deformed. Then, the fabric sheet is allowed to drape on the sphere. Figure 17c–f depict the top, side and isometric views of the final draped shape.

123

Closure Grid-based methods for fabric drape simulation often resort to simple particle mechanics and formulate the elastic energy in terms of the change in the inter-particle distance and the change in the angle between the straight lines joining adjacent particles. The complexity of these irrational quantities increases with the order of the differentiation which is required in both the optimization and the Newton–Raphson solution procedures. In this paper, gridlines along warps and wefts are setup on a fabric sheet. The gridlines will enable the fabric sheet to be discretized into integration zones. The stretching strain is quantified by the normal Green strain. By virtue of a corotational frame, small strain and small curvature assumptions, the curvature of a gridline segment bounded by three nodes can be approximated by a linear combination of the nodal displacement vectors.

Int J Mech Mater Des (2005) 2:225–243

241

(b)

(a) Y

X

Z

Y

X

Z

(d)

(c)

X X

Y

Y Z

Z

(e) (f) Fig. 16 Two 92 · 180 cm fabric sheets jointed by a 2 cm wide seam along the weft direction

With the small strain assumption, the shear strain can be approximated as the dot product of the two deformed gridline segments that define the integration zone. Consequently, all the elastic energies commonly considered in drape simulations are simple polynomials of the nodal displacement in

the present model. Fabric drape problems with extensive wrinkles are then examined. The presence of sharp fold, seam and cut in the undeformed fabric sheet as well as fabric-to-solid contact are considered. The predicted appearances are pleasant and conform to real life observations.

123

242

Int J Mech Mater Des (2005) 2:225–243

Fig. 17 A sharp fold along the weft direction in the middle of a 180 · 180 cm fabric sheet Acknowledgment The financial support of the Research Grant Council of Hong Kong(project number : HKU 7204/ 03E) is gratefully acknowledged.

References ABAQUS.: ABAQUS Theory Manual Version 6.1. Hibbitt, Karlsson & Sorensen Inc. Rhode Island (2000)

123

Ascough, J., Bez, H.E., Bricis, A.M.: A simple finite element model for cloth drape simulation. Int. J. Cloth. Sci. Tech. 8, 59–74 (1996) Breen, D.E., House, D.H., Getto, P.H.: A particle-based model of woven cloth. Visual Comput. 8, 264–277 (1992) Chen, B., Govindaraj, M.: A physical based model of fabric drape using flexible shell theory. Textile Res. J. 65, 324–330 (1995)

Int J Mech Mater Des (2005) 2:225–243 Choi, K.J., Ko, H.S.: Stable but responsive cloth. SIGGRAPH 2002 Conference Proceedings, Annual Conference Series: ACM Press/ACM SIGGRAPH: pp. 604–611 (2002) Dai, X.Q., Furukawa, T., Mitsui, S., Takatera, M., Shimizu, Y.: Drape formation based on geometric constraints and its application to skirt modeling. Int. J. Cloth. Sci. Tech. 13, 23–37 (2001) Eberhardt, B., Weber, A., Strasser, W.: A fast, flexible particle-system model for cloth draping. IEEE Comput. Graphics Appl. 16, 51–59 (1996) Eischen, J.W., Deng, S., Clapp, T.G.: Finite element modeling and control of flexible fabric parts. IEEE Comput. Graphics Appl. 16, 71–80 (1996) Etzmuss, O., Keckeisen, M., Strasser, W.: A fast finite element solution for cloth modelling. Proceedings of the 11th Pacific Conference on Computer Graphics and Applications, pp 244–251. IEEE (2003) Gan, L., Ly, N.G., Steven, G.P.: A study of fabric deformation using nonlinear finite elements. Textile Res. J. 65, 660–668 (1995) Gong, D.X., Hinds, B.K., McCartney, J.: Progress towards effective garment CAD. Int. J. Cloth. Sci. Tech. 13, 12–22 (2001) Hinds, B.K., McCartney, J., Woods, G.: 3D CAD for garment design. Int. J. Cloth. Sci. Tech. 4, 6–14 (1991) House, D.H., DeVaul, R.W., Breen, D.E.: Towards simulating cloth dynamics using interacting particles. Int. J. Cloth. Sci. Tech. 8, 75–94 (1996) Kang, T.J., Yu. W.R.: Drape simulation of woven fabric by using the finite-element method. J. Textile Inst. 86, 635–648 (1995)

243 Lin, J.Z., Tang, M., Dong, J.X.: Cloth simulation based on local adaptive subdivision and merging. The 8th International Conference on Computer Supported Cooperative Work in Design Proceedings, pp. 718–726, IEEE (2003) Meyer, M., Debunne, G., Desbrun, M., Barr, A.H.: Interactive animation of cloth-like objects in virtual reality. J. Visual. Comput. Animation 12, 1–12 (2001) Rodel, H., Ulbricht, V., Krzywinski, S., Schenk, A., Fischer. P.: Simulation of drape behavior of fabrics. Int. J. Cloth. Sci. Tech. 10, 201–208 (1998) Stylios, G., Wan, T.R., Powell, N.J.: Modeling the dynamic drape of garments on synthetic humans in a virtual fashion show. Int. J. Cloth. Sci. Tech. 8, 95–112 (1996) Sze, K.Y.: A corotational beam model for fabric drape simulation. Book of Abstract, keynote paper, 4th International Conference for Mechanics & Materials in Design, Nagoya, Japan, 4–7 June 2002 (2002) Sze, K.Y., Liu, X.H.: A corotation grid-based model for fabric drapes. Int. J. Numer. Meth. Eng. 57, 1503–1521 (2003) Tan, S.T., Wong, T.N., Zhao, Y.F., Chen, W.J.: A constrained finite element method for modeling cloth deformation. Visual Comput. 15, 90–99 (1999) Teng, J.G., Chen, S.F., Hu, J.L.: A finite-volume method for deformation analysis of woven fabrics. Int. J. Numer. Meth. Eng. 46, 2061–2098 (1999) Villard, J., Borouchaki, H.: Adaptive meshing for cloth animation. Eng. Comput. 20, 333–341 (2005) Zhang, D., Yuen, M.M.F.: Cloth simulation using multilevel meshes. Comput. Graphics 25, 383–389 (2001)

123