Selected Methods of Obtaining Adequate Forms of Deformation of Supercritical Thin-Walled Flat Plates

The study considers the numerical methods of obtaining the correct forms of deformation of thin-walled plates subjected to pure shear, operating under the conditions of critical loads. Selected methods were presented and subjected to comparative tests, allowing for the obtaining of a state of large deformations of a perfectly flat system, loaded in its plane. These methods are based on the application of various forms of initial imperfections into the system, either geometrically or physically. A new solution was proposed in this area, differing in its essence from most commonly used solutions.


INTRODUCTION
Among the many types of load-bearing structures used in modern technology, a significant role is played by thin-walled structures, commonly used when it is necessary to meet rigorous mass criteria. In most cases, their correct operation is expected in the load ranges that do not cause the effects of stability loss, and the buckled elements of the structure are treated as damaged. Only some types of thin-walled structures, such as aircraft structures, are an exception.
Most of the typical air load-bearing structures are based on the semi-monocoque model, the unquestionable advantage of which is a relatively low weight in relation to the possible load capacity of the system. A feature that distinguishes the aviation semi-monocoque structures from most other structures is the limited acceptability of post-critical deformations. This phenomenon is permissible only in relation to the fragments of coverings, limited by skeleton elements (e.g. ribs and stringers), and therefore it must have a local character. In the case of a post-critical deformation of any of the remaining structure elements, the element is considered to be damaged. The elastic nature of the phenomenon is also an obligatory requirement.
Controlled loss of stability of shell fragments in aircraft structures is allowed primarily for metal structures, although in some cases it is also applicable for composite structures (e.g. loadbearing structures of the Boeing 787 and Airbus 350 aircrafts). While several types of aeronautical load bearing components are characterized by a pronounced shell curvature (e.g., fuselage skin segments), many others are flat systems. An example of a design scheme with completely planar components is a wing spar with a wall forming a thin-walled structure subjected to shear (Fig. 1).
Similar solutions are used in various parts of the airframe structure, e.g. in the case of a wing spar, but also in the case of the parts of the fuselage or tail elements. However, the fragments of the fuselage and wing coverings are also characterized by relatively small curves in the case of large aircraft. With a sufficiently dense arrangement of the skeleton elements, from the computational point of view, they can be treated as flat. In all of the shell systems mentioned, the phenomenon of loss of stability due to shear stresses can Selected Methods of Obtaining Adequate Forms of Deformation of Supercritical Thin-Walled Flat Plates occur. In many cases, this phenomenon is also acceptable, either in the range of permissible loads or in the range between their permissible and destructive values. In most cases, it is necessary to perform advanced supercritical deformation analysis in the design processes. In the flat segments of the covering, the tangential stress distributions caused by shear have the character of a drawing field, which under real system conditions are the cause of geometric imperfections, as a result of which the supercritical deformations appear without the occurrence of distinct bifurcations and in the case of rectangular shells, they usually have the form of oblique folds (Fig. 2).
An important problem related to the numerical mapping of this type of deformation in the case of non-linear numerical analyses using the finite element method is the natural lack of tendency to the formation of the imperfections mentioned, resulting from the idealized geometric form of the model. Therefore, it is necessary to take them into account in the form of an appropriate extortion, of a geometric or physical nature.
In this study, the advantages and disadvantages of selected methods used in this field were analysed, and a solution that would allow obtaining correct results of numerical analyses with relatively easy adjustments to the basic system, based on a perfectly flat model, was proposed.

PURPOSE AND SCOPE OF THE STUDY
The subject of the considerations was a thinwalled square plate, subjected to pure shear. With its geometry and boundary conditions, it corresponded to the system, which is the subject of the planned experimental studies, using a model made of polycarbonate (Fig. 3). The edges of the plate are stiffened with steel cladding, articulated in the corners. The fastening system consisted of two supports, located in two opposite corners (Fig. 4). The first one (point A) deprived the joint of translational degrees of freedom and rotation about the axis u (Fig. 4). The second (point B) blocked the displacements corresponding to the axis directions The value of the load applied in the lower corner of the system was 500 N. The numerical representation of the diagram presented above was based on the use of surface elements in relation to the tested thin-walled structure and eight-node three-dimensional elements in relation to the stiffening edges. The choice of the latter solution was dictated by the desire to maximally simplify the model geometry by not using contact functions. In the adopted geometry variant, the articulated joints were modelled by using the appropriate shape of the stiffening corners, using the properties of a three-dimensional element, characterized by the lack of active rotational degrees of freedom (Fig. 5).
The models use a finite element mesh of a regular nature and a density ensuring a satisfactory convergence of the solution in the linear range. The mesh consisted of 21,830 three-dimensional elements and 21,138 shell elements with linear in the anticipated experimental system For numerical analyses in all cases, the MSC Patran/MSC MARC software was used, which is characterized by an extensive implementation of numerical methods that enable the mapping shape functions (Fig. 6). The aim of the considerations was to compare the effectiveness of the selected methods of obtaining the correct form of supercritical deformations of the analysed object under the conditions of non-linear numerical analysis. For this purpose, several selected variants of the system based on the assumptions regarding geometry, boundary conditions and loading presented above were analysed.

NUMERICAL ANALYSES
Numerical mapping of advanced deformation states in the case of flat systems, loaded in their plane, is related to the problem resulting from the idealized geometric form of the model and the direction of load application. In the case of real systems, the presence of natural imperfections is the factor that initiates the appearance of displacements in the direction normal to the plane of the object. In the case of an idealized model, the result of the nonlinear numerical analysis is an idealized distribution of displacements occurring only in the model plane (Fig.7).
Therefore, in order to obtain correct results of the supercritical deformation analysis, it is necessary to use one of the methods of applying imperfections of a geometric or physical nature.
For some types of commercial software based on the finite element method, it is possible to apply the procedures to implement the results of a linearized stability analysis as data used by the nonlinear analysis. In such cases, the first part of the numerical problem involves the application of the method based on the Euler stability criterion, also known as the method of adjacent states. As a result, it leads to solving the problem on eigenvalues: where:K 0 is the stiffness matrix resulting from the material properties, determined for the initial state, and λ • K 1 stiffness matrix resulting from the current geometry of the system, depending on the state (load) control parameter l. The essential nonlinear analysis, in turn, consists in determining the course of the equilibrium path of the system, which in the general case is a hypersurface in n-dimensional state space, where n is the number of degrees of freedom of the system. In each subsequent step of the analysis, related to the load increment, the discrete system must satisfy the matrix equation of residual forces: r(u, λ) = 0 (2) where: u is a state vector containing the displacement components of structure nodes corresponding to its current geometric configuration, λ is a control parameter corresponding to the current load level, and r is vector residual containing the unbalanced force components related to the current state of deformation of the system. A feature of all nonlinear procedures is the presence of an incremental phase. For each successive increment, at the transition from state n to the state n + 1, the values not specified are the increments: In order to determine them, an additional equation is formulated, called the increment control equation, or constraint equation: c(Δu n , Δλ n ) = 0 (4) constituting a condition defined by the user, resulting from the adopted correction strategy.
In the discussed case, due to the procedural limitations resulting from the applied software type, on the basis of the lowest form of the model buckling determined in a linearized manner, another geometric model was prepared, containing the resulting imperfections. In the next, non-linear step, using the Newton-Raphson method and the correction based on the state control, the supercritical deformation distributions were determined, corresponding to the assumed load value (Fig. 8).
Although the linearized pre-buckling analysis in many cases is an unreliable method and does not allow for the correct determination of the buckling mode, in the analysed case, due to the nature of the stress distribution corresponding to the drawing field, the obtained results can be considered as reliable in terms of quality. A precise quantitative analysis would only be possible with the results of the experiment. Therefore, recognizing the correctness of the presented solution, it should be emphasized that, as it was mentioned, in the case of a number of types of commercial software, the method used is sometimes labourintensive and does not allow reproducing the deformation states corresponding to the load ranges slightly exceeding the critical values.
The aim of these considerations is to present an alternative method that allows initiating a geometric imperfection in a relatively simple way, which in turn makes it possible to numerically map supercritical deformations, already the load level corresponding to the critical value. The method is based on the application of a load in the form of concentrated forces applied at the corners of the model edge stiffening system, in the direction normal to its median plane (Fig. 9).
It causes a slight deflection of the system in the u-v plane, as a result of which the numerical model loses its perfectly flat geometric character. This effect is a desired imperfection  in the case of analyses, causing the initiation of local bifurcations in the state hyperspace. The forces in the drawing marked with the symbols P n by assumption, have the values significantly lower than the basic load applied in the model plane. An important factor determining the correctness of the solution seems to be the application of their correct values. In the next step of numerical analyses, a model was used without preliminary geometric imperfections, additionally loaded with the forces corresponding to the diagram presented above. As expected, the form of supercritical deformations turned out to be largely dependent on the value of the exciting forces (Fig. 10).
As a result of a series of analyses, during which successive increases in the exciting load were applied, a diagram was drawn up illustrating the relationship between the ratio of the force value P n to the value of the load in the plate plane and the value Δλ corresponding to the difference between the absolute value of the maximum displacement λ in the direction perpendicular to the model plane and the value considered to coincide with that obtained via the previously described method using linearized buckling analysis λ 0 , referred to λ 0 (Fig. 11). for too low value of the exciting force, b) solution similar to that obtained in the previous analysis step Fig. 11. Graph of the dependence of Δλ and the force ratio P n to the load value The above-mentioned dependence proves that there is a very narrow range of values of the force causing imperfection, at which it is possible to obtain a correct solution. The rapid increase in the value of Δλ corresponds to a very rapid increase in the height of the fold on the main diagonal of the model with a relatively small increase in imperfection. In turn, a sharp decrease in said value with a further increase in the value of P n corresponds to an increase in the height of the side folds with a simultaneous decrease in the fold height on the main diagonal (Fig. 11).

CONCLUSIONS
The presented method of initiating the process of loss of stability of the analysed thin-walled structure by introducing transverse loads of small values, resulting in a geometric imperfection, seems to be effective, despite its simplicity. Its additional advantage is the potential possibility of tracing the initial phase of supercritical deformation, which does not seem possible in the case of specifying the buckling form resulting from the linearized pre-buckling analysis.
The disadvantage of the method in relation to the aforementioned comparative procedure is a very small range of values of transverse forces causing the imperfection, allowing for obtaining a correct solution. In the case when the distribution of displacements corresponding to supercritical deformations is determined by experiment or calculations carried out by another method, recognized as reliable, the selection of the values of transverse forces may be carried out with the method of successive approximations until a satisfactory solution is reached. It cannot be ruled out that the relationship between the values of the exciting forces and the load in the plane of the model is fixed, either generally or for a specific type of material. However this thesis must be proven based on careful analyses of a number of models, supported by a model experiment.