Coupling of zones with different resolution capabilities in dynamic finite element models of woven composites/ Skirtingos skiriamosios gebos zonu sujungimas austu kompozitu dinaminuose baigtiniu elementu modeliuose.

Rimavicius, V. ; Calneryte, D. ; Barauskas, R. 等

1. Introduction

Finite element modelling of material and structure physical behaviour have been already an integral part in the development of new products, and as well a part of design processes for civil engineering, mechanical and mechatronics engineering, munitions industry, and in many other engineering areas. Numerical models of biomechanical and biomedical engineering systems are applied in the analysis of implant mechanics, of dental prosthesis strategies, in the development of protective clothes and lightweight armour, planning of surgery and hyper-thermal treatment procedures of tumours, etc. Computer simulation results often are highly adequate to the reality and reliable in such a high degree that nowadays many natural experiments may be replaced by calculations.

The finite element computational technology enables to investigate highly refined structural models, as well as, complex models of material response in order to represent the behaviour of real systems accurately and adequately, including their essentially non-linear behaviour, characteristics of inner structure, peculiarities of internal texture etc., [1,2]. Though powerful computers and perfect computational software are widely available nowadays, real needs of engineering investigation steadily require more computational resources than is available at the moment. The search of reasonable compromise between the accuracy of the model and its dimensionality is termed as model reduction.

Reduced models are created essentially for every theoretical and computational investigation. It is always a choice of a researcher to decide what minuteness of structure details or of its integral structural elements should be represented by the model. Very often the choice is based on engineering intuition and proper understanding of the physical essence of the behaviour of the investigated system. However, nowadays the availability of computational technologies implies new regular approaches of treatment of the same phenomena at different time and length scales (multiscale modelling). The same object can be described by using different mathematical abstractions. The change of the model resolution property is usually carried out by generating elements based on mathematical equations presented in different forms and by introducing necessary interactions rather than by coarsening or refining the finite element mesh.

The model resolution requirement has been recognized as one of the main challenges in non-linear computational mechanics almost two decades ago [3]. The problem remains actual nowadays as well, though a lot of research have been carried out under multiscale analysis headline. The reason is obviously a tremendous variety of manifestations of the multiscale. The term model making has been proposed in [4] and demonstrated that the building a proper computational model could be understood as a certain kind of engineering activity. Also term model engineering appeared [5]. Reduction concepts and techniques of mathematical models, which appear in very different application areas, such as dynamic analysis in multiple time scales, functional dependence extraction, chemical kinetics, population dynamics, etc., have been discussed in [6-7].

In this work, along the physically-based analogies used for dynamic interaction model reduction two formal approaches have been investigated. A model earlier developed by us in LS-DYNA finite element system has been employed as a reference model. The internal structure of the reduced models synthesized in this work is essentially simpler as it presented the woven structure by means of the girder, the interactions between nodal points of which were expressed by using elastic, viscous and friction link finite elements. Contact interactions of the surface of the bullet against the girder needed to be determined during each time integration step. The mathematical expression of mutual adequacy of the two models was based on the comparison of the response of the reduced model against the response of the reference model by means of appropriately constructed penalty functions.

2. Finite element model of a woven structure

As an example of a complex dynamic interaction problem, a model of impact of a homogeneous deformable bullet against the multilayer woven textile structure has been analysed (Fig. 1).

Each fabric layer is made of multifilament yarns woven together. We used the "mezzo-mechanical" approach, in which a yarn has been considered as primary component comprising the fabric. The properties of an individual yarn are defined empirically by making several assumptions. It s commonly accepted that the cross-section of a free yarn not interwoven into a textiles is close to a circle made of cross-sections of a certain number of individual filaments.

In a fabric structure the yarns are being compressed because of forces acting in overlapping areas. As a result, the geometry of the cross-section of each yarn is being changed depending upon the constitution of the yarn, its density, its type and technological parameters of the weave. After examination of the cross-section of the yarn extracted from the fabrics we assume it to be close to the combination of two circular segments. On the base of the measured cross-section thicknesses at nodal points, shell element structure has been employed in order to represent each yarn of the weave. As shell thicknesses are prescribed and may not change during the contact interaction, under certain interaction conditions this may affect the adequacy of the model.

[FIGURE 1 OMITTED]

The technique of obtaining the model of a woven yarn patch in more detail has been described earlier in our earlier publication [8]. By using this technique the weave is obtained by solving the contact interaction equilibrium problem and is able to present the warp and weft yarns in the weave. The Ox directed yarns (the warp) are elastically deformed by moving their nodes in the direction Oz perpendicular to the plane of the layer in order to situate the Oy directed yarns (the weft) in-between the warps. After activation of the contact search algorithm the yarns of the model come to an elastic equilibrium. As a result, a woven structure is obtained in a similar way as a real fabric is produced. After equilibrium is obtained, the elastic stresses and strains are being fully or partially removed in order to imitate the relaxation of stresses in multi-filament yarns. Moreover, the amount of residual stress may be controlled when building the model. In this way we may create a desired value of the warp tension, which is an important characteristic of a weave.

3. Multi-scale model based on the adjustment of wave propagation properties

The zone of single layer of a fabric taking place in the contact interaction for a 9 mm bullet comprises only about 20-60 mm in diameter. However, in order to represent properly the dynamics of the interaction process much larger pieces of a fabric should be modelled. Unfortunately, the weave step being about 1 mm, such "fully woven" models are prohibitive because of their huge dimensionality. In order to obtain models of smaller dimensions the macro-mechanical approach is being used. The zones of the fabric remote from the interaction zone are modelled by means of orthotropic shell (membrane) elements. At the zones distant form the interaction point that can be much bigger than the elements, which present the mesh of interacting yarns. Sudden transitions in mesh density of finite element models are implemented by using tied interfaces, Fig. 2.

[FIGURE 2 OMITTED]

A uniform membrane-type surrounding of the woven patch cannot be expected to present an identical dynamic behaviour as the woven yarn structure. However, satisfactory approximation is possible provided that the appropriate selection of geometrical and physical properties of the membrane zone is made. Two models of the fabric layer with different dimensions of the woven patch are built, Fig. 3. The model containing a larger patch has been considered as a reference model (in our study it contained 120x120 yarns), and the smaller one was 20x20 yarns; we take it at least twice greater than the diameter of the zone at which the failure of yarns may take place.

The membrane is presented by 1-integration point shell elements exhibiting no bending stiffness. The membrane thickness may be assumed as known, however, the values of mass density p, Young's modulus E and shear modulus G of the orthotropic material have to be selected in order to obtain the dynamic response of both models as close as possible to each other.

The overall approach is similar to the convergence investigation of models of different refinement. The "exact" numerical solution may be regarded as known as we are able to make a stand-alone numerical experiment by employing the fully woven model. Further a series of models with different sizes of the woven patch are investigated and the physical constants of the surrounding orthotropic membrane established in order to ensure a convergent behaviour in the sense that the solutions obtained by using every model are close to the exact one.

In order to estimate the wave propagation speed as the wave propagates in between two nodes of the structure we fix time moments [t.sub.1], [t.sub.2] when the displacements of the nodes reach a pre-selected level, see Fig. 4, a and b.

[FIGURE 3 OMITTED]

For instance, the average speed of the wave propagating along direction x may be estimated as v = [increment of x]/[t.sub.2] - [t.sub.1] where [increment of x] = [x.sub.2] - [x.sub.1] distance between the nodes. The discontinuities in the wave propagation speed relationships are eliminated at each nodal point by taking the average of the propagation speeds in two neighbouring inter-nodal segments of the wave propagation line.

The results presented in Fig. 3 are obtained by using one of "successful" sets of membrane material constants as

Surface density of the membrane/surface density of the woven patch = 0.8;

Young's modulus of the membrane/Young's modulus of the yarn = 0 3.

Shear modulus of the membrane/Young's modulus of the membrane = 0.0042.

Rather small values of Young's modulus of the equivalent membrane enables to present the de-crimping driven deformation of the fabric which can be expected to prevail in zones remote from the point of impact. Very small values of shear modulus are obtained because of small shear stiffness of the woven structure at small strains. However, the value of shear modulus had to be selected with care as it makes a significant influence basically upon the shape of the front of the transverse wave, as well as, upon its amplitude. It should be noticed that the wave propagation speed relationships of the reference model and 10x 10 woven patch model have significant (up to ~20%) differences.

The resultant displacements and wave fronts obtained by using 10x 10 model are close to the results obtained by using the reference model as can be observed from the results presented in Fig. 4.

[FIGURE 4 OMITTED]

4. Obtaining reduced finite element models of the woven patch

As a demonstration of the technique a simple problem has been formulated, in which the rectangular girder model composed of interconnected rods has been reduced to equivalent membrane model. Consider structures of a girder and equivalent membrane [9]. Both structures are characterized by their physical and geometrical properties. Suppose that the membrane is described by the following set of physical parameters: Young's module [E.sub.m], Poisson's ratio [u.sub.m], shear module [G.sub.m], thickness [s.sub.m], material density [[rho].sub.m] and number of elements Nm along one side of the structure. The girder is described by using Young's module [E.sub.g], material density [[rho].sub.g], rod width [b.sub.g], rod thickness [h.sub.g] and the number of cells [N.sub.g] along the side of the girder. For the first approach we select uniform finite element models of the same size for both the girder and membrane and same positions of nodes. Both models are subjected to the same loads F and the same boundary conditions. A subset of reference nodes is selected and marked by circles in Fig. 5. Reference node displacements p of equivalent membrane and q of the girder are obtained by solving the finite element analysis problem.

[FIGURE 5 OMITTED]

Displacements of the girder depend on the girder parameters vector x = {[b.sub.g], [N.sub.g]}, therefore q = q(x). Displacements p = p(y) of membrane nodes depend on the vector y = {[G.sub.m], [s.sub.m], [[rho].sub.m] } of membrane physical parameters. Values of membrane elements physical parameters y are to be found, which ensure minimum mismatch of displacements between reference nodes of equivalent membrane and girder finite element models. For evaluation of mismatch we define the error function T(y, x) = T(p(y), q(x)). Since p(y) and q(x) are nonlinear functions, T(y, x) also is non-linear function. We formulate the error function constrained minimization problem as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

where [[M.sup.q]], [[M.sup.p]] , [[C.sup.q]] , [[C.sup.p]] , [[K.sup.q]] , [[K.sup.p]] are the mass, damping and stiffness matrices of the girder and membrane structures correspondingly; {q}, {p} are the nodal displacement vectors; {[F.sup.q]} , {[F.sup.p]} are the nodal force vectors; [x.sub.nx1] is the state parameters vector; [y.sub.nx1] is the optimization parameters vector, [a.sub.i] [less than or equal to] [y.sub.i] [less than or equal to] [b.sub.i], i = [bar.1,n]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

where [p.sub.i]. = {[P.sub.ix],[P.sub.iy]} and [q.sub.i] = {[q.sub.ix],[q.sub.iy]} , i = [bar.1,k], are displacement vectors of reference nodes, t - time, 0 [less than or equal to] t [less than or equal to] T. As sample physical values in this investigation [E.sub.m] = [E.sub.g] = 1.5x[10.sup.11] Pa, [[upsilon].sup.m] = 0, [b.sub.g] = [h.sub.g], [p.sub.g] = 7800 [k.sub.g]/[m.sub.3], are assumed. The girder rods were assumed to satisfy the length-thickness relation as:

1/20 [less than or equal to] [b.sub.g]/L[less than or equal to]1/8 [??] L/20[N.sub.g] [less than or equal to] [b.sub.g] [less than or equal to] L/8[N.sub.g] (3)

Error function minimization was carried out in the static and dynamic modes for different sets of girder physical parameters. The sets of membrane physical parameters have been found, which ensured the minimum mismatch between the displacements of both models.

As an alternative approach the neural network (ANN) method has been employed. Two-layer neural network has been used, Fig. 6, a.

Back-propagation learning algorithm was used for training. In order to avoid the overtraining of the network the Early Stopping algorithm was employed.

The neural network parameters n1 and f1 have been sought, which produced the smallest value of the error function. A single hidden-layer was used with neural number in the range 3 [less than or equal to] [n.sub.1] [less than or equal to] 17 and transfer function [f.sub.1] = {"tansig", "logsig"}. 500 experiments were performed. The smallest error of test vectors was [MSE.sub.m]in [approximately equal to] 0.03333, where [n.sub.1] = 11 and [f.sub.1] = "tansig", Fig. 6, b. By using this technique functional dependences between given vectors of physical parameters of the girder x = {[b.sub.g], [N.sub.g]} and equivalent membrane y = {[G.sub.m], [s.sub.m], [[rho].sub.m]} were obtained [9], which are further discussed in Section 5.

As the further model reduction step the rectangular girder and equivalent membrane zones within the combined model have been considered, where mesh densities of the girder and of the membrane were different, Fig. 5. The geometrical parameters of the equivalent membrane of the combined FE model presented in Fig. 5, b have been employed, which ensured that the central zone of the model presented by the dense mesh of interconnected rectangular rods of dimension [L.sub.m] x [L.sub.m] responded like the central zone of the girder in Fig. 5, a subjected to the same mechanical forces [10].

[FIGURE 6 OMITTED]

Optimization problem (1) has been solved by regarding all the nodes within the [L.sub.m] x [L.sub.m] central zone of the girder as reference nodes. Functional dependences were established between pairs of physical parameter sets obtained after solution of the optimization task. The regression and, as an alternative, the neural network approaches were applied. Numerical experiments were carried out and the comparison between results obtained by using the two approaches was performed [10], see the discussion in Section 5.

[FIGURE 7 OMITTED]

The structural behaviour of the reduced models in the dynamic was analysed by introducing time as a continuous variable t [member of][0;T]. The optimization was performed

by using the target function in form (2), which was integrated in time over time interval [0; T]. Nodal displacements were stored at all time steps [t.sub.k].

The duration of the investigated time interval was selected as T [approximately equal to] L [([E.sub.g]/[[rho].sub.g]).sup.-1/2], which is necessary for the longitudinal elastic wave to travel distance L equal to the side length of the model. The sine-pulse shaped time law of the excitation force F had the period [T.sub.F] [approximately equal to] T/2 . Forces at in dividual nodes have been applied as F = 2 [F.sub.max] sin ([pi]t)/2 N-1

Functional dependence between physical parameters of the girder and equivalent membrane x = {[b.sub.g],[N.sub.m],[N.sub.g]} and y ={[G.sub.m],[s.sub.m],[[rho].sub.m],[DELTA]} was obtained, where [DELTA] is the estimation of the mismatch between the girder and the equivalent membrane reference nodes displacements).

The final model reduction task was formulated in order to simplify a complex finite element computational model, which imitated the ballistic contact interaction of a bullet against the textile structure implemented in LSDYNA finite element software. The numerical solution obtained in LS-DYNA environment was regarded as the reference solution. The reference computational model was developed and verified earlier on the base of real experiments [8]. Here we simplified the task of the reduced model development as the failure of the yarns was not considered during the bullet against woven structure interaction.

The computational model was reduced by replacing the woven structure by the girder model, Fig. 7. The relative movement of the rods during the deformation of the structure was represented by interconnecting the rods by means of visco-elastic links at intersection points, Fig. 8.

[FIGURE 8 OMITTED]

The interconnection links of the rods at their intersection points can be presented as structural link elements described as follows: where [k.sub.x],[k.sub.y],[k.sub.z] are stiffness elements along the directions Ox,Oy,Oz correspondingly, [c.sub.x],[c.sub.y],[c.sub.z] are damping (viscous friction) elements along the directions Ox, Oy, Oz correspondingly.

As in the previous problems, the error function was based on the nodal displacements of reference and reduced FE models integrated over time interval [0; T ] and by considering all the nodes of equivalent girder rods intersections as reference nodes. The stiffness and damping constants [k.sub.x], [k.sub.y], [k.sub.z] and [c.sub.x], [c.sub.y], [c.sub.z] of the link elements were regarded as optimization variables of the error function minimization task.

5. Computed results

The rod-based girder structure 11 x 11 in the quarter symmetry model has been created. It was surrounded by the equivalent membrane zone presented by ANSYS elements SHELL63. Mesh roughness of the membrane zone was increased with the distance from the girder zone, which enabled to reduce the total number of elements. The unilateral constraint elements CONTA175 and TARGE170 were employed for presenting the contact interaction of the woven structure against the bullet. The surface of the bullet was assumed as absolutely rigid and was implemented by SHELL63 elements.

Primarily the stiffness coefficients of the link elements at connection points of girder rods were determined by using the trial-error approach as [k.sub.x] = [k.sub.y] = 3.8661x[10.sup.6] N/mm in x and y directions and [k.sub.z], = 5.2448 x [10.sup.3] N/mm in z direction. The error function between responses of the reference and reduced FE models was based on the cumulative differences of displacements of appropriate nodes, Fig. 2, a and Fig. 9, a.

Displacements and errors of displacements of reference nodes of the reference and reduced FE models along the xy line are compared in Fig. 10.

In Fig. 10 the numbers are accompanied by the symbols "k" and "e", where "k" indicates that the node belongs to the reduced FE model and "e" indicates that the node belongs to the reference model.

In Fig. 10, a the reference nodes displacements in both models can be observed as they change in time. The maximum error is approximately 47% at nodes 6393, 6831 and 6719 of the reduced model, Fig. 10, b.

[FIGURE 9 OMITTED]

[FIGURE 10 OMITTED]

The error function minimization problem has been solved in order to adjust the values of the parameters of the reduced model. The MATLAB function fmincon has been applied in order to solve the error function minimization problem (1), (2).

Coefficients [k.sub.x], [k.sub.y], [k.sub.z] and [c.sub.x], [c.sub.y], [c.sub.z] of the link elements were regarded as optimization variables. It was assumed that [k.sub.x] = [k.sub.y] because of the quarter symmetry of the reduced FE model. The constraints on the values of the stiffness coefficients were imposed as:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

The solution demonstrated that the coefficient values [k.sub.x] = [k.sub.y] = 3.8661x[10.sup.3] N/mm, [k.sub.z] = 5.2448x[10.sup.6] N/mm provide acceptable results.

Maximum error between reference nodes of the reference and reduced FE models was 14% after the adjustment of the stiffness coefficient values of link elements, Fig. 11.

[FIGURE 11 OMITTED]

6. Conclusions

The issues of the reduced model development for essentially non-linear structural interactions have been discussed basing on the example of the finite element analysis of the interaction process of the woven structure against a rigid projectile. The reduced model was synthesized by using highly refined finite element reference model. The basic tasks and algorithms were programmed in MATLAB. The finite element software ANSYS and LSDYNA were used for obtaining the structural responses at given sets of model parameters.

The size of the woven structure reference model was reduced to reasonable dimensions by using a mezzo-mechanical approach, in which the multifilament yarns comprising the woven fabric were presented as a structure of narrow bands of the prescribed cross-section and by presenting the zones of the fabric distant from the point of impact by the orthotropic membrane model. The physically-based approach of combining together the models of woven structure and the equivalent membrane structure was developed by employing the coincidence of in-plane and out-of-plane wave propagation velocities in both zones as a prerogative for the equivalency between the two differently scaled models.

The development of the reduced models included approaches based on the cumulative error function minimization in both static and dynamic modes of the non-linear structural analysis. The reduced finite element model has been developed, which tended to replace a complex interaction model of mechanical contact among woven yarns by the linear or close-to-linear model. The girder structure interconnected by appropriate structural link elements has been employed. The mechanical contacts had to be analyzed only in the rigid projectile against the yarns interactions thus enabling essential reduction of the structural complexity of the model. The regression and artificial network techniques were involved for the synthesis of the reduced models, which were designed to represent properly the structural behaviour provided no yarns failure takes place.

As demonstrated by the presented examples, the presented approach was able to present satisfactory results both in static and dynamic modes of structural analysis. Nevertheless, in some cases the direct application of the approach may encounter difficulties as the error function may appear as multi-extreme and difficult to minimize globally. On the other hand, the approach inevitably starts with certain physically-based knowledge about the essential features of the analyzed system. Very probably the model reduction approach can never be "purely formal". Though based on particular examples, the presented analysis creates some methodical prerogatives for the reduction of models of other systems as well and in this way promotes to the field of research of multi-scale models.

10.5755/j01.mech.19.34654

Received October 26, 2011 Accepted April 08, 2013

Acknowledgements

The research has been sponsored by Lithuanian Science Council, agreement number MIP-31/2011.

References

[1.] Zienkiewicz, O.C. 2000. Achievements and some unsolved problems of the finite element method, International Journal for Numerical Methods in Engineering, 47(1): 9-28. http://dx.doi.org/10.1002/(SICI)1097-0207(20000110/30)47:1/3<9::AID-NME793>3.0.CO;2-P.

[2.] Zienkiewitcz, O.C.; Taylor, R.L. 2000. The Finite Element Method fifth edition. Publishhed by Butter worth-Heinemann, ISBN 0750650494.

[3.] Belytschko, T. 1996. On Difficulty Levels in Non Linear Finite Element Analysis of Solids, IACM Expressions, Bulletin of the International Association for Computational Mechanics, 2: 6-8.

[4.] Peierls, R. 1980. Model-making in physics, Contemporary Physics, 21: 3-17. http://dx.doi.org/10.1080/00107518008210938.

[5.] Gorban, A.N.; Kazantzis, N.; Kevrekidis, I.G.; Ottinger, H.C.; Theodoropoulos, C. 2006. Model Reduction and Coarse-Graining Approaches for Multiscale Phenomena. Springer, Berlin. http://dx.doi.org/10.1007/3-540-35888-9.

[6.] Barth, T. J. et al. (Eds). 2010. Coping with Complexity: Model Reduction and Data Analysis, Springer, p.348.

[7.] Fedaravicius, A.; Kilikevicius, S. et al. 2011. Research of the mine imitator interaction with deformable soil and its practical realisation, Mechanika 17(6): 615617.

[8.] Barauskas, R., Abraitiene, A. 2007. Computational analysis of impact of a bullet against the multilayer fabrics in LSDYNA, International Journal of Impact Engineering 34(7): 1286-1305. http://dx.doi.org/10.1016/jijimpeng.2006.06.002.

[9.] Barauskas, R.; Rimavicius, V. 2008. Reduction of finite element models by parameter identification, Information Technology and Control 37(3): 211-219.

[10.] Barauskas, R.; Rimavicius, V. 2010. Coupling of zones with different resolution capabilities in finite element models of uniform structures, Information technology and centrol 39(1): 7-17.

V. Rimavicius *, D. Calneryte **, R. Barauskas ***

* Kaunas University of Technology, Studentu str.50-407, 51368 Kaunas, Lithuania, E-mail: vidmantas. rimavicius@gmail.com

** Kaunas University of Technology, Studentu str.50-407, 51368 Kaunas, Lithuania, E-mail: d.calneryte@stud.ktu.lt

*** Center for Physical Sciences and Technology, Demokratu 53, 48485 Kaunas, Lithuania, E-mail: rimantas.barauskas@ktu.lt