Numerical Reproducing of a Bifurcation in the Stress Distribution Obtaining Process in Post-Critical Deformation States of Aircraft Load-Bearing Structures

1. Introduction

Modern aviation structures are characterised by widespread application of thin-shell load-bearing systems. The strict requirements with regard to the levels of transferred loads and the need to minimise a structure mass often become causes for accepting physical phenomena that in case of other structures are considered as inadmissible. An example of such a phenomenon is the loss of stability of shells that are parts of load-bearing structures, within the range of admissible loads.

Thus, an important stage in design work on an aircraft load-bearing structure is to determine stress distribution in the post-critical deformation state. One of the tools used to achieve this aim is nonlinear finite elements method analysis. The assessment of the reliability of the results thus obtained is based on the solution uniqueness rule, according to which a specific deformation form can correspond to one and only one stress state. In order to apply this rule it is required to obtain numerical model’s displacements distribution fully corresponding to actual deformations of the analysed structure.

An element deciding about a structure’s deformation state is the effect of a rapid change of the structure’s shape occurring when the critical load levels are crossed. From the numerical point of view, this phenomenon is interpreted as a change of the relation between state parameters corresponding to particular degrees of freedom of the system and the control parameter related to the load. This relation, defined as the equilibrium path, in case of an occurrence of mentioned phenomenon, has an alternative character, defined as bifurcation. Therefore, the fact of taking a new deformation form by the structure corresponds to a sudden change to the alternative branch of the equilibrium path [1-4].

Therefore, a prerequisite condition for obtaining a proper form of the numerical model deformation is to retain the conformity between numerical bifurcations and bifurcations in the actual structure. In order to determine such conformity it is required to verify the results obtained by an appropriate model experiment or by using the data obtained during the tests of the actual object. It is often troublesome to obtain reliable results of nonlinear numerical analyses and it requires an appropriate choice of numerical methods dependent upon the type of the analysed structure and precise determination of parameters controlling the course of procedures.

Due to the number of state parameters, the full equilibrium path should be interpreted as hyper-surface in state hyperspace, satisfying the matrix equation for residual forces:

$$r(u,\Lambda )=0,$$E1

where u is the state vector containing structure nodes’ displacement components corresponding to current geometrical configuration, is a matrix composed of control parameters corresponding to current load state, and r is the residual vector containing uncompensated components of forces related to current system deformation state. The set of control parameters may be expressed by a single parameter that is a function of the load. Equation (1) takes then the following form:

$$r(u,\lambda )=0,$$E2

called a monoparametric equation of residual forces.

The prediction-correction methods of determining the consecutive points of the equilibrium path used in modern programs contain also a correction phase based on the satisfaction of an additional equation by the system, called an increment control equation or constraints equation:

$$c(\Delta u_n,\Delta {\lambda }_n)=0,$$E3

where the increments:

$$\Delta u_n=u_{n+1}-u_{n}and\Delta {\lambda }_n={\lambda }_{n+1}-{\lambda }_n$$E4

correspond to the transition from n-th state to n+1-th state.

The graphic interpretation of the increment control equation is presented in Figure 1.

In order to find out whether there is full conformity between the character of actual deformations and their numerical representation it would be required to compare the combinations of the relevant state parameters in all the phases of the course of the phenomenon considered herein. Because of the complication of such a comparative system, the deformation processes are represented in practice by applying substitute characteristics called representative equilibrium paths. They define the relations between a control parameter related to load and a selected, characteristic geometric value related to a structure’s deformation, an increment of which corresponds to a change of the value of all or some state parameters [5-9].

Figure 1.

Graphic presentation of various correction strategies for the representative system of one degree of freedom

In case of a large number of state parameters it is not possible at all to represent the character of bifurcation by applying a representative equilibrium path. Sometimes, changes of state parameters resulting from local bifurcation may show the lack of perceptible influence on the representative value, which results in non-occurrence of any characteristic points on the representative path. In general, however, these changes cause a temporary drop in the control parameter value (Figure 2).

Figure 2.

Bifurcation points on a representative equilibrium path (u - representative geometric value,λ – control parameter related to load).

So both the experiment itself and nonlinear numerical analyses may result only in a representative equilibrium path. In that case, the problem of the numerical representation of bifurcation comes down to the preservation of conformity of the representative equilibrium path obtained by a numerical method with the one obtained experimentally, where a sine qua non for the application of the solution uniqueness rule is to recognise the similarity of the post-critical deformation forms of the experimental and numerical models as sufficient [10,11].

An additional problem occurred during the experimental determination of the equilibrium path, resulting from the lack of abilities of recording the said temporary, little drops in load, arising from local bifurcations, causing changes of the values of some state parameters. In the majority of experiments, the load of the tested model is achieved by force control, e.g. using a gravitational system, or displacement, by means of various types of load-applying devices (Figure 3).

Figure 3.

A stand for testing thin-shell structures subject to torsion: left – a version with a loading system controlling the displacement (turnbuckle), right – a version with a system controlling by force (gravitational).

However, even in case of devices with high level of technical advancement, in general it is not possible to register precisely short-lasting force changes, occurring from the beginning of a bifurcation phenomenon to the moment of reaching the consecutive deformation form by the model. Therefore, the representative equilibrium path obtained as a result of the experiment is of smooth characteristics, and its formation is based on measuring points corresponding to the consecutive deformation states determined.

In case of nonlinear numerical analysis in the finite element approach, the accuracy of the obtaining of the representative equilibrium path may be much more accurate. The existing commercial programs usually offer the results of all the increment steps, followed during the calculation process, and thus they also allow observing slight fluctuations of the control parameter. The only limitation here is exclusively the value of the incremental step itself. In spite of this, due to the lack of possibilities of relating the results obtained to the relevant detailed changes of the experimental characteristics, it seems appropriate to determine the numerical representative equilibrium path of the same level of simplification as in the case of the experiment.

2. Analyses of example structures

The comparative analysis of such representative equilibrium paths is not, however, a method that allows a complete enough verification of the reliability of the results of numerical calculations. An example of a problem in which the calculated results have been deemed incorrect despite the seeming full conformity of the representative equilibrium paths is a thin-shell open cylindrical structure with edges strengthened by stringers, working in the conditions of constrained torsion (Figure 4). This type of systems is quite often used in aviation structures. They form areas of cockpits and large cut-outs, e.g. in cargo airplanes and they are usually adjacent to much stiffer fragments of the structures [12,13].

Figure 4.

A schematic view of the tested structure and comparison of the deformation forms obtained experimentally and numerically

The area adjacent directly to the closing frame turned out to be crucial in the problem under consideration. The stringer strengthening the edge of the structure was buckled, and the experimental model sustained plastic deformation. The relation between the total torsion angle of the examined structure and the torque moment constituting the load was adopted as the representative equilibrium path. In spite of the seeming conformity of the obtained characteristics, a different character of deformation was observed in the critical area (Figure 3). The divergence seems to result here from the symmetric character the bifurcation phenomenon initiating deformation. Different deformations correspond to two possible variants of the actual equilibrium path, the differentiation of which is not possible in case of using a simplified representative path. So, in spite of the conformity of the characteristics, the effective stress distributions obtained numerically cannot be considered as reliable.

The need for an analysis of stress distributions in case of structures similar to the one presented above occurs quite rarely due to the commonly adopted principle, pursuant to which beam structures after buckling are considered as damaged [14]. The considerations relating to shells used in the aviation industry, e.g. semi-monocoque structure elements are of much greater practical significance.

Examples of such a system are open cylindrical shells, which were subjected to a cycle of tests, during which it was assumed that stringers are characterised by a sufficient margin of stiffness and they do not lose stability (Figure 5).

Figure 5.

Geometry of the model

So the models made of polycarbonate were strengthened with longitudinal members with large, rectangular cross-section and relatively high values of geometric moments of inertia. Each of the examined systems was subjected to constrained torsion, using a test stand presented in Figure 3. The tests aimed at developing the methodology of determining stress distributions in a structure’s shell in post-critical deformation states.

In the conditions of torsion the stress state created in the shell of such system, may be interpreted as an incomplete tension field. As a result, even if there are no geometrical imperfections, the shell loses its stability. On the other hand, the post-critical deformation increment causes significant stress redistribution. The experiments repeated a number of times showed that the final form of post-critical deformations of such systems, occurring at sufficiently high load values, is always the same in spite of the alternative character of the course of the structure state changing process (Figure 6).

Figure 6.

A post-critical deformation form of a two-segment structure at maximum load (the total torsional torque equals 180Nm)

The fact that local bifurcations following increases in loads occur with some scatter of locations and stress levels makes the nonlinear numerical analysis particularly troublesome in this case. It is practically impossible to develop a FEM model allowing to reproduce accurately the entire process of the structure’s state changes, using commercial software, due to the nature of the functioning of algorithms for choosing the variants of the equilibrium path at the bifurcation points and the impossibility of the user’s interference in the form of those algorithms. In this situation it seems appropriate to focus only on obtaining a numerical solution consistent with the experiment results at assumed load values.

The selection of an appropriate combination of numerical methods and parameters controlling the course of the analysis seems particularly vital in this case, likewise the proper representation of the model’s stiffness. Even small mistakes in this respect result in the occurrence of incorrect forms of deformation (Figure 7).

Figure 7.

Incorrect forms of shell stability loss, showing a member’s buckling not revealed in the experiment

It should be emphasised that it seems very risky to rely in the design process on the results of nonlinear numerical analysis of similar structures without appropriate verification in an experiment, if only a relatively cheap model experiment. In practice, multiple repetition of the analysis and systematic comparison of its results with the results of the experiment are required to obtain correct results of the numerical representation of a structure’s state in the conditions of post-critical loads (Figure 8).

Figure 8.

Accepted as satisfactory form of post-critical deformation obtained as a result of nonlinear analysis (left) and a deformation form obtained as a result of experimental tests at identical load (right)

Examples of aviation load-bearing structures, being much less troublesome from the point of view of nonlinear numerical analyses, working in the conditions of post-critical deformations are semi-monocoque box structures with flat walls. This type of solutions is often applied in modern military aircrafts. They usually constitute compact load-bearing structures of fuselages and internal wing areas, most frequently ogival or delta-shaped.

A representative element of this type of structures that was subjected to a detailed experimental and numerical analysis was a multi-segment multi-member thin-shell structure, working in the conditions of dominant torsion (Figure 9).

Several variant of the presented structure were examined, including the ones with openings of various shapes, corresponding to all types of functional and service cut-outs.

Striving to preserve the local character of post-critical deformations of the segments of skin, similar structures possess relatively thin shells at considerable stiffness of the framing. So, a characteristic feature of the structure’s shell is a low level of critical load. Even at little levels of shearing stress, the waving of the skin occurs within particular segments. At the same time, as proved by the experiments carried out, to reach the global loss of stability corresponding to the damage of the structure it is required to increase the load value considerably. Thus, in spite of the natural susceptibility to a local loss of skin stability, the structure is characterised by very favourable properties from the point of view of their application in aviation.

Figure 9.

The examined representative case of a box structure

The process of the appearance of consecutive bifurcations is in this case of established and repeatable character. This results from the mere mechanism of tension fields formation in particular segments of the skin. Additionally, the character of the stress distribution in each of the structure’s segments is similar. All these factors result in the fact that in contrast to the structure discussed earlier, nonlinear numerical analysis of semi-monocoque box structures is possible when using alternatively various numerical methods, and the successful course of the analysis does not require the application of correction strategies based on arc length control methods (Figure 1).

The numerical analysis of the examined structure allowed to obtain results of high degree of conformity with the results of the experiment (Figure 10).

The conformity between the courses of the representative equilibrium paths did not raise any objections either, and their form itself constitutes the confirmation of a gentle nature of the phenomenon (Figure 11).

While analysing the constructional solutions used in aviation load-bearing structures, the need should be emphasised for drawing attention to the area of working shells, in which local cracks occur due to cyclic loads. Such damages are usually dimensionless in nature, with edges touching each other. In case the skin is stretched along the normal to the direction of the crack, in the vicinity of the front of the crack there occur both tensile stresses and compressive stresses (crosswise), which creates conditions for the occurrence of a local buckling. In these conditions, in the weakened zone there occurs a phenomenon of a loss of shell stability within the area of a crack, referred to as wrinkling [11].

Figure 10.

Comparison of the deformation form in a structure with a circular opening obtained during the experiment (left) and as a result of MES nonlinear analysis (right)

Figure 11.

Comparison of the representative equilibrium paths – a structure with a circular opening

During the experimental and numerical analysis of the phenomenon the model shown in Figure 12 was used.

In order to examine the influence of a crack length on the nature of the phenomenon, a number of the models versions were analysed. The quantitative and qualitative character of deformations in the vicinity of the crack was determined by means of the Shadow moiré method (Figure 13).

Although during the experiment with many repetitions of load application, in the form of a force stretching the examined plate along the normal to the direction of the crack edges (Figure 11), the post-critical deformation was repeatable in nature, it turned out troublesome to obtain correct results of nonlinear numerical analyses. The deformation phenomenon is relatively gentle in nature, which results from the susceptibility of the structure to the slightest geometric imperfections near the crack.

Figure 12.

The model for experimental examination of the wrinkling phenomenon

Figure 13.

The Shadow moiré patterns for three crack lengths: 30mm, 50mm, 70mm

According to the expectations, an attempt to carry out a nonlinear numerical FEM analysis using an ideally flat plate model resulted in the settlement of the state of equilibrium of the structure which did not show any displacements in the perpendicular direction to the plate surface. In such a case, bifurcation point does not appear at all on the equilibrium path in the numerical space of the state.

The nature of the analysis undergoes a complete qualitative change in case of taking imperfection into consideration, either in the form of initial displacements of crack edges, or in the form of preload residual loads. The degree of the representation of the nature of the phenomenon depends mostly on the type of the finite element. As results from the numerical test performed, very few elements offered by commercial software manufacturers possess a sufficiently appropriate mathematical model of displacements to represent in full the phenomenon of wrinkling. They are usually elements of rich, nonlinear shape functions, the application of which in complex nonlinear tasks causes in turn the lack of effectiveness of algorithms for choosing the proper branches of the equilibrium path. They may result in e.g. the occurrence of the so-called “ping-pong effect”, consisting in reversing the course of the analysis and its further advancement within the range of negative values of the control parameter.

From the numerical point of view, the mere complexity of nonlinear procedures forces the application of a possibly uncomplicated displacement model of the used finite elements. In case of elements of linear shape functions, recommended by the manufacturers of MSC MARC software for nonlinear analyses, the correct selection of the type and size of the initial imperfection turned out to be crucial for forcing the appropriate course of the phenomenon.

It turned out the most effective to apply continuous cross-bending load of constant intensity along both edges of a crack, possessing the value by several orders lower in relation to the force loading the plate. The post-critical deformation distributions obtained as a result of numerical analyses were compared with the experiment both in qualitative and quantitative terms, which was possible owing to the moiré patterns.

After numerous repetitions of the analyses and the selection of the most appropriate set of numerical methods, a satisfactory displacement distribution was obtained (Figure 14).

Figure 14.

Comparison of deflection distributions in the post-critical state for a model with a crack 30 mm long

Owing to the examination of models with cracks of various lengths, it became possible to determine the nature of relation between the crack length and the resultant value of cross-bending load, constituting a factor causing initial imperfection (Figure 15).

Figure 15.

The relation of the values of total cross-bending load in the crack length function

The research results of various load-bearing structures confirm that the difficulty related to carrying out an appropriate nonlinear numerical analysis results from the nature of bifurcation. If the change in a structure’s form is gentle in nature and it occurs in a small area, then the bifurcations related to it occur gradually, in relatively small subsets of state parameters. The numerical simulation of the process is then easy to perform and it may take place when using prediction-correction methods with simple correction based on state control. But if the deformation occurs in a larger area, and the change of the form is violent in nature, then the bifurcation corresponds to the simultaneous change of a great number of state parameters, and the determination through a numerical procedure of their appropriate combination, corresponding to the new state of static equilibrium, may be hindered or even impossible. In such a case it is necessary to apply matchings of prediction methods with correction strategies based on arc length control methods, such as the Riks correction or the Crisfield hyperspherical correction [15] (Figure 9).

Bifurcation changes of the forms of load-bearing structures, containing shells of considerable curvature, occur more violently if there is a higher relation of the square of the smaller of dimensions of the shell segment area limited by the adjacent members’ frames to the value of the local radius of its curvature [16]. Thus, semi-monocoque structures of relatively low number of the framing elements are especially troublesome in nonlinear numerical representations.

An example of such a structure is a closed cylindrical shell presented in Figure 16.

Figure 16.

A schematic view of a complete cylindrical shell reinforced by four members (a) and a schematic view of a structure including dimensions (b)

The structure’s framing consists of a minimum number of crosswise elements, i.e. two closing frames and four longitudinal members. The type of the structure itself corresponds to solutions commonly used in the aviation technology, e.g. the construction of a fuselage of an aircraft. It should be emphasised, however, the model subjected to examinations constitutes a special instance of a structure of purposefully minimised number of longitudinal members. The actual solutions are usually based on much more extended framings. The structure described corresponds to an isolated phase of a wider cycle of examinations aiming at determining direct dependences between the number of framing elements and post-critical deformation distributions.

The examined structure was subjected to constrained torsion using a modified version of the stand presented in Figure 2.

According to the expectations, post-critical deformations occurred in a violent way. Due to the gravitational way of load application, the measurement of the relation between torsion angle and the torque moment, assumed as the representative equilibrium path, corresponded to the steady states (Figure 21).

Using this mode of taking measurements, the representative characteristics does not reflect bifurcation points in an overt way, but attention should be drawn to the occurrence of its horizontal section. It corresponds to this phase of the experiment in which a sudden change occurred in the structure state with the simultaneous constant load level.

With regard to the symmetry, the deformed structure possessed four characteristic grooves in all the shell segments (Figure 17). During the experiment the surface geometry was registered using the projection moiré method. Atos scanner manufactured by a German company, GOM Optical Measuring Techniques was used as a registering device.

Figure 17.

The advanced post-critical deformation of the examined structure (left) and the distribution of contour lines representing the size of the deformation, made using the projection moiré method (right)

The problem discussed belongs to one of the most troublesome from the point of view of a FEM nonlinear numerical simulation. A number of tests performed using the MSC MARC software revealed the lack of effectiveness of its procedures in case of this problem, with regard to determining the appropriate post-buckling state of a structure. The algorithms used in those procedures are characterised by inability to represent the symmetry of the phenomenon. With the idealised geometric form of the model, the obtaining of the new form of the structure after crossing the critical load value occurs only in one of the segments, in spite of the apparently correct, symmetrical initiation of stability loss. This proves the faults in the algorithms for choosing the appropriate variants of the equilibrium path in case of the appearance of changes in the state parameters combination in several of their independent subsets.

The situation was improved when shell imperfections were implemented, by applying normal forces to the skin, in the central points of particular skin (Figure 18).

Figure 18.

A geometrical model of a structure made in MSC PATRAN environment with boundary conditions and loading

However, even in the case of applying this type of forcing a form change, it was very difficult to obtain results that would fully correspond to the experimental results. Assuming the use of skin elements with linear shape functions, the appropriate density of the mesh turned out to be the key factor, but its excessive density caused incorrect forms of post-critical deformations (Figure 19).

Figure 19.

The incorrect form of deformation, obtained in case of too many elements

The better result, in case of application of beam elements as a representation of stringers, was obtained with the use of a relatively low density of mesh. This proves the rightness of the thesis, proved a number of times in many studies, pursuant to which the decrease in the general number of degrees of freedom, corresponding to the number of state parameters, in case of nonlinear procedures used in the available commercial programs, often brings benefits that considerably exceed the deficiencies of a mathematical description resulting from the decrease in the number of elements.

The best result was obtained only after the fundamental change of the concept of FEM model, when the different kind of finite elements was applied as a representation of stringers (thick shell element was used instead of a recommended beam element). However this solution, from the point of view of mathematical description is much less correct, it turned out much more effective in case of relatively low values of the total torsion angle of the structure.

The results of analysis of this FEM model version, obtained using secant prediction method and strain correction strategy (Figure1), are presented in Figure 20.

Figure 20.

The deformation distribution (left) and reduced stress distribution acc. to Huber-Mises hypothesis (right) for 100% of the maximum load (stringers modelled with thick shell bilinear element)

The strain-correction strategy turned out most effective in case of significant, violent change of the form of deformations, when the representative equilibrium path contains relatively long “horizontal section”.

The relation between the representative equilibrium paths is presented in Figure 21.

Figure 21.

The presentation of the representative equilibrium paths

3. Conclusion

The presented examples of load-bearing structures represent only some of those many used in the modern aviation technology. But the criterion applied while selecting them as objects of experimental and numerical analyses was its representativeness for the most commonly met elements of constructions, in case of which the occurrence of a local stability loss is acceptable in the conditions of service load.

The fundamental conclusion that can be drawn from the presented results of the research is the absolute need for using experimental verifications with regard to FEM nonlinear numerical analyses of this type of structures. The more so that even in the cases in which the correctness of the results obtained seems unquestionable, they may be in fact burdened with errors resulting from the very limited reliability of the numerical procedures used in commercial programs.

Based on the nonlinear numerical analyses, related to the presented structures, frequently repeated many times, a general recommendation may also be formulated for the maximum possible limitation of the size of a task. Striving for increasing the accuracy of the calculations by increasing the density of finite elements mesh, applied successfully in linear analyses, may turn out ineffective in case of a nonlinear analysis and may lead to incorrect results or the lack of convergence of calculations.

The numerical representation of bifurcation, by virtue of the mere idea of the discrete representation of continuous systems, must be simplified in case of the finite elements method. In such a situation, based on the quoted examples, the need must be emphasised for obtaining the indispensable convergence of the experimental and obtained numerically relations between a selected geometric parameter characterising the essence of a structure’s deformation and a selected value relating to the load, recognised as representative equilibrium paths. This convergence, in combination with the accepted as sufficient similarity of post-critical deformation forms, constitutes the grounds for accepting the reliability of stress distributions determined by means of numerical methods.