Coupling of zones with different resolution capabilities in dynamic finite element models of woven composites

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


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 com-putational 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.

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.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.

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 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 120120 yarns), and the smaller one was 2020 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 , 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 1 , t 2 when the displacements of the nodes reach a pre-selected level, see Fig. 4 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 1010 woven patch model have significant (up to ~20%) differences.
The resultant displacements and wave fronts obtained by using 1010 model are close to the results obtained by using the reference model as can be observed from the results presented in Fig. 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 struc-tures 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 m , Poisson's ratio υ m , shear module G m , thickness s m , material density  m and number of elements N m along one side of the structure.The girder is described by using Young's module E g , material density  g , rod width b g , rod thickness h g and the number of cells N 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.Displacements of the girder depend on the girder parameters vector } , { where

As sample physical values in this investigation
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 n 1 and f 1 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 ≤ n 1 ≤ 17 and transfer function   L  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.
The structural behaviour of the reduced models in the dynamic was analysed by introducing time as a continuous variable   where Δ 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 LS-DYNA finite element software.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.The interconnection links of the rods at their intersection points can be presented as structural link elements described as follows: where ,,

Computed results
The rod-based girder structure 11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 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. .k .
The solution demonstrated that the coefficient

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 synthe-sized by using highly refined finite element reference model.The basic tasks and algorithms were programmed in MATLAB.The finite element software ANSYS and LS-DYNA 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 mezzomechanical 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.

COUPLING OF ZONES WITH DIFFERENT RESOLUTION CAPABILITIES IN DYNAMIC FINITE ELEMENT MODELS OF WOVEM COMPOSITES
S u m m a r y An approach for creating a reduced cost-effective dynamic finite element computational model of the missile impact upon the woven textile target is presented.The model of a woven patch developed in LS-DYNA finite element environment, the internal structure of which was based on the mathematical representation of all contact interactions among the yarns of the woven structure has been employed as a reference model.The equivalent membrane model served for the computationally cheaper representation of the farther zones of the woven structure.They farther zones surround the patch of the woven structure, which is in immediate contact with the missile.The synthesis of the membrane model is based on the adjustment of in-plane and out-of-plane waves propagation velocities in both woven and membrane zones.The reduced model of the woven structure is based on the girder structure, the links between nodal points of which were expressed by using elastic and viscous elements.Only contact interactions of the surface of the missile against the girder needed to be determined during each time integration step.The task of the synthesis of the reduced model has been presented in the form of the optimization problem, where the cumulative displacement error has been used as the target function.The physical meaning of the target function is the mismatch between the displacements of the two models end.Alternatively, the artificial neural network approach has been employed for obtaining the physical parameters of the reduced model.

Fig. 1
Fig. 1 Quarter-symmetry model of the lead bullet-fabric package interaction

Fig. 2
Fig. 2 Fragments of woven fabric quarter symmetry models with indicated wave propagation directions and numbers of nodes for comparison of dynamic properties: a)woven patch 1010 yarns; b)fragment of the woven patch 6060 yarns (reference model)

Fig. 3
Fig. 3 Time laws of longitudinal (a) and transverse (b) displacements of nodes and determining impact wave propagation speeds along direction x in the vicinity of the tied interface For instance, the average speed of the wave propagating along direction x may be estimated as

Fig. 4
Fig. 4 The shape of the transversal wave front and zdisplacement contour plot obtained by using reference model (a) and 1010 combined model (b) (at time instant 610 -5 s)

Fig. 5 FE
Fig. 5 FE models of girder (a) and rectangular membrane with girder zone and reference node (b)

F
 K are the mass, damping and stiffness matrices of the girder and membrane structures correspondingly;   q ,   p are the nodal displacement vectors;   q

Fig. 6
Fig. 6 Two-layer neural network (a); average square error dependence on network parameters (b) Optimization problem (1) has been solved by regarding all the nodes within the m m L 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 k .The duration of the investigated time interval 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 2 dependence between physical parameters of the girder and equivalent membrane

Fig. 7 Fig. 8
Fig. 7 Reference model (a) and reduced FE model (b) 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 ,, of the link elements were regarded as optimization variables of the error function minimization task.
in z direction.The error function

Fig. 9 Fig. 10
Fig. 9 Reduced FE model with indicated reference nodes (a) and displacements of the reduced FE model at time moment 5 Fig. 11 Displacements (a) of model reference nodes on the xy diagonal line (Fig. 9, a) and their errors (b)