Nonlinear bifurcation analysis of stiffener profiles via deflation techniques Thin-Walled Structures

When loading experiments are repeated on different samples, qualitatively different results can occur. This is due to factors such as geometric imperfections, load asymmetries, unevenly stressed regions or uneven material distributions created by manufacturing processes. This fact makes designing robust thin-walled structures difficult. One numerical strategy for exploring these different possible responses is to impose various initial imperfections on the geometry before loading, leading to different final solutions. However, this strategy is tedious, error-prone and gives an incomplete picture of the possible buckled configurations of the system. The present study demonstrates how a deflation strategy can be applied to obtain multiple solutions for the more robust design of thin-walled structures under displacement controlled uniaxial compression. We first demonstrate that distinct initial imperfections trigger different sequences of instability events in the Saint Venant – Kirchhoff hyperelastic model. We then employ deflation to investigate multiple bifurcation paths without the use of initial imperfections. A key advantage of this approach is that it can capture disconnected branches that cannot easily be discovered by conventional arc-length continuation and branch switching algorithms. Numerical experiments are given for three types of aircraft stiffener profiles. Our proposed technique is shown to be a powerful tool for exploring multiple disconnected bifurcation paths without requiring detailed insight for designing initial imperfections. We hypothesise that this technique will be very useful in the design process of robust thin-walled structures that must consider a variety of bifurcation paths.


Introduction
The investigation of nonlinear bifurcation in thin-walled structures, especially for imperfection-sensitive shells where multiple bifurcation paths are possible [1], makes the bifurcation analysis and design particularly challenging. These issues have been classically addressed by means of artificially triggering one specific bifurcation path using geometric imperfection. The use of linear buckling analysis is not suitable because of the nonlinear degradation of membrane stiffness of the shells when imperfections are present [2]. Discrepancy between the linear theory and the experimental observation is well known from the classical buckling literature [3][4][5][6][7][8][9][10] and also verified in dynamic analysis of cylindrical shells [11].
Castro et al. [12] compared five different imperfection patterns that are commonly used, the most common one being the use of linear buckling modes as imperfections (abbreviated LBMIs) in the nonlinear buckling case. This is frequently applied in civil engineering to analyse tanks, silos, and cooling towers under various load conditions. Essentially, the structure is perturbed using a small multiple of an eigenmode obtained with linear buckling analysis. Studies on imperfection sensitivity using LBMIs can be found in the works by Yamada and Croll [13,14]. Sosa et al. [15] used LBMIs to create the initial prescribed displacement field required in the reduced energy method. Schmidt [16] and Winterstetter and Schmidt [17] applied nonlinear bifurcation analysis for steel towers under wind loading. In the field of space launcher structures, Hilburger et al. [18] compares the effect of LBMIs on cylinders with imperfections coming from manufacturing signatures. LBMIs are also known for their use in triggering nonlinear buckling bifurcation paths in many other works [6,[19][20][21][22][23][24][25][26][27][28][29]. An advantage of LBMIs over measured imperfections is that the imperfection pattern is easily obtained and included in a finite element model [12,28].
Another method that gained recent attention to induce geometric imperfection in thin-walled shells consists of simply using perturbation loads [12,25,27,28,, with the aim of creating single buckle-like imperfection patterns. Using LBMIs or perturbation loads, the amplitude and shape of the imperfection pattern is chosen by intuition, and the nonlinear analysis performed for each perturbation will follow only one bifurcation path. In the present study, we demonstrate how the deflated continuation algorithm [53,54] can greatly increase the robustness of nonlinear bifurcation analysis in the sense that no choices have to be made regarding the perturbation to follow different bifurcation paths, nor about their corresponding amplitudes. This algorithm is different from the conventional combination of arc-length continuation and branch switching algorithms [54]. The arc-length continuation method has been applied to evaluate nonlinear bifurcation dynamics in thin-walled structures, such as stiffeners with low torsion stiffness [55,56], plates [57,58] and shells [11,26,57,[59][60][61].
The combination of arc-length continuation and branch switching can robustly trace out the branches connected to a specified initial solution. However, they do not identify branches that are disconnected from the given data without further intervention. This disconnection frequently arises in applications due to geometry asymmetry of the design; continuation in the symmetry of the geometry is not generally feasible. In Ref. [62], Cox et al. introduce a new structural optimisation technique, called modal nudging, by perturbing the perfect baseline geometry to nudge the systems onto higher load-carrying paths. This allows the exploration of disconnected bifurcation paths and shows great potential in helping to design structures with improved load-carrying capacity, compliance and stability. However, information about buckling modes is required for this technique. In contrast, the deflation technique [54] to be applied in this work does not require any a priori information about buckling modes while still enabling the discovery of disconnected paths.
Three practical aircraft stiffeners are considered: the L-shaped asymmetric type, the L-shaped symmetric type and the Z-shaped asymmetric type, as detailed in Fig. 1. It is important to emphasise that, in the present study, the nonlinear post-buckled bifurcation paths are investigated for the stiffener profiles in a non-assembled configuration, i.e. not as part of a stiffened panel. In real applications, the stiffeners would be assembled in a panel, where post-buckling configurations of such stiffeners can usually be achieved before the ultimate loads supported by the structure.
All stiffeners are modelled with Saint Venant-Kirchhoff hyperelastic materials and used to test the algorithm under displacement-controlled uniaxial compression. All implementations are performed in Firedrake [63] and rely on the PETSC [64] and the Defcon library [65].
The remainder of this paper is organised as follows. We briefly review the Saint Venant-Kirchhoff hyperelastic model in Section 2. We then illustrate how the imposition of geometric imperfections can trigger different bifurcation paths for the aircraft stiffeners in Section 3. In order to capture branches that may be disconnected due to the asymmetry of the stiffener designs, we review the deflation technique and demonstrate its application in bifurcation analysis of the three stiffener profiles in Section 4. Finally, conclusions are drawn in Section 5.

Saint Venant-Kirchhoff hyperelastic model
Consider a three-dimensional body occupying a reference configuration B 0 with Lipschitz continuous boundary Γ subject to certain loads to the body, thus leading to a deformed configuration B . In this model, we characterise the deformation by the displacement u : B 0 → R 3 and define the deformation gradient tensor by where I is the identity second-order tensor. For Saint Venant-Kirchhoff hyperelastic materials, the constitutive equation (i.e., the stress-strain relation) can be written as (2.1) with λ and μ being Lam� e parameters, σ the second Piola-Kirchhoff stress tensor and E the Lagrangian strain tensor given by In order to have a well-posed problem, additional boundary conditions are needed. We divide the boundary Γ into two disjoint parts: with the Dirichlet boundary Γ D and the traction boundary Γ N . On the top and bottom boundary faces with ε being a parameter that will be continued in the deflated continuation algorithm [54]. Note that the z-direction corresponds to the longitudinal (axial) direction of stiffeners in this work. For simplicity, we assume that the aircraft stiffeners are homogeneous, isotropic and frame-indifferent [66].
Remark 2.2. The displacement-controlled boundary condition on Γ D guarantees that all degrees of freedom at Γ D have the same displacement condition. Imposing a traction on the bottom boundary may not yield an even displacement distribution.

Remark 2.3.
Note that ε in the boundary data u 0 corresponds to the axial displacement applied to the stiffeners.
On Γ N , we have no traction, i.e., where the traction tðuÞ is defined by tðuÞ ¼ PðuÞn; with P ¼ σ F denoting the first Piola-Kirchhoff stress tensor and n the outward normal to the boundary surface. In addition, the uniaxial applied force f ext can be calculated from the second Piola-Kirchhoff stress tensor σ via Here, the negative sign is added to represent the positive compressive force. Then, it is known that the average stress over the bottom face is computed by with jΓ bottom j denoting the measure of the bottom face. The boundary value problem considered in this work is with b being the body force vector. In the implementation, we let B 0 be the aircraft stiffener and ignore the gravitational body force, i. e, b ¼ 0, as it is negligible compared with the compressive force that we impose. Denote the admissible function space of the displacement by The weak form of (2.7) can be derived as: The Dirichlet boundary condition u ¼ u 0 will be enforced weakly later using Nitsche's method.
Remark 2.4. The W 1;4 -regularity is needed to make the weak form (2.8) well-defined. Indeed, by direct computations, we can obtain If u;v 2 W 1;p , then PðuÞ 2 L p=3 and rv 2 L p . Thus, PðuÞ : rv is in L p=4 . This requires p ¼ 4 at least for (2.8) to be well-defined. Moreover, by the Sobolev embedding theorem [66], Theorem6.1-3], the W 1;4 -regularity guarantees that pointwise evaluation is well-defined.

Enforcement of the essential boundary condition
The traction-free boundary condition (2.4) is naturally enforced in (2.8) via the divergence theorem; it remains to enforce the essential boundary condition (2.3). Throughout this work, we will follow Nitsche's method [67] to weakly impose the Dirichlet boundary condition u ¼ u 0 on Γ D . To this end, we add the following two terms to the weak form (2.8). Here, γ > 0 is a large penalty parameter, necessary for numerical stability [68]. Note that the second term in (2.9) arises from integration by parts of the divergence of the first Piola-Kirchhoff stress tensor P in (2.7).
Remark 2.5. The term R ΓD tðuÞ⋅v ds in (2.9) is well-defined. Indeed, this can be seen from the inequality [69,70].
with c > 0 a mesh-dependent constant. Consequently, we summarise the final variational problem used in this work as follows: find u 2 V such that (2.10) Essentially, (2.10) is a consistent formulation as the additional penalty term is zero for an exact solution.
Remark 2.6. Here, we use the non-symmetric version of Nitsche's method [71] for ease of the stability analysis (see Appendix A).
Furthermore, we can see that the variational problem (2.10) is nonlinear due to the presence of the nonlinear stress-strain relation (2.1). Hence, the classical Newton method is applied and the r-th Newton iteration takes the form of with the update δu 2 V to the current approximation u r . Here, the first Gâteaux derivative is given by where with the fourth-order elasticity tensor C in the form of Here, δ ij is the Kronecker delta defined by and A ij denotes the ði; jÞ-th entry of the second-order tensor A.
In the implementation, we take the parameters E ¼ 69 GPa, ν ¼ 0:334 for aluminium stiffeners [72,73] and choose γ ¼ 10 15 based on unreported preliminary experiments. The choice of γ is related to the stability analysis of the Newton system (2.11), discussed in Appendix A. Further discussions of the Saint Venant-Kirchhoff model and other models can be found in Refs. [66,74].

Geometric imperfections inducing different bifurcation paths
Real manufactured profiles will show different imperfection patterns that may exhibit various mechanical responses. The present discussion aims to show the limitation of current nonlinear algorithms that are capable of only capturing one bifurcation path depending on the initial imperfection pattern.
The following analyses were performed using NX Nastran and conventional Newton-Raphson iterative schemes, a tetrahedral mesh with mesh size of 1.3 mm and continuous piecewise quadratic Lagrange elements for the displacement. The use of solid elements is preferred due to the high-fidelity discretisation of the cross-section and the direct compatibility of the mesh with the technique to be used in the sequel. Fig. 2 shows the effect of a perturbation load of 20N on the bifurcation path of the L-shaped asymmetric profile, compressed up to 0.6 mm. When multiple bifurcation paths exist, the instabilities that happen first decide the fate of the succeeding ones, strongly affecting the load at which the remaining component will fail. Conventional solvers are only able to evaluate one of such sequences of events. Fig. 3 shows the axial compression of a Z-shaped asymmetric profile up to 2 mm. Note that in the case without the perturbation load the right-hand side flange undergoes a local buckling, whereas the application of a perturbation load of 20N changes the critical buckling flange to be the left-hand side one, as indicated in Fig. 3. Results show that even a small initial disturbance may decide which sequence of bifurcation paths will be followed by conventional nonlinear analysis solvers. For the L-shaped profile with a bulb and the Z-shaped profile herein evaluated, a perturbation load of only 20N is enough to change the bifurcation path.
These study cases with L-and Z-shaped profiles just illustrate what generally happens in optimised thin-walled structures. Other larger examples are wing and fuselage structures of real aircraft, which are highly optimised so that the buckling of different components usually happen at similar loads [75][76][77]. In these designs, small imperfections due to variations in manufacturing parameters can induce different bifurcation paths and hence different buckling loads [78]. For instance, if one decides to perform ultimate failure tests in different aircraft wings, the region where it fails can be different for each case, because the slightest disturbance might be sufficient to trigger a different bifurcation path. Designers can control the sequence of instability events by selectively increasing the margins of safety in some regions, distancing the bifurcation events from one another. However, increasing the margins of safety will result in heavier designs that should be avoided. Therefore, more robust numerical techniques that enable the investigation of multiple bifurcation paths are required. This will significantly increase the chance of predicting real test results, even in structures with very close bifurcation paths, which is generally the case in thin-walled shells. We will return to this issue in the next section by introducing the deflation technique.

Deflated continuation methods
As discussed in the previous section, more robust methods investigating multiple bifurcation paths should be considered. We notice that due to the presence of non-linearity in the Saint Venant-Kirchhoff model (2.10) and possibly geometric imperfections in our concerned aircraft stiffeners, the variational problem (2.10) can permit multiple equilibrium states. In this section, we will briefly review the deflation technique and use the deflated continuation algorithm proposed in Ref. [54] to exploit the buckling profiles for three practical types of stiffeners (see Fig. 1).

Deflation
Consider a parameter-dependent nonlinear problem f ðu; λÞ ¼ 0 for u 2 U and λ 2 R; (4.1) where U is an admissible space for u and λ is the parameter. For our purposes, we assume that this problem permits multiple solutions for some values of λ, which we wish to find. Its bifurcation diagram will then visualise how solutions change as the parameter λ varies over the range ½λ min ; λ max �.
A classical strategy used in solving (4.1) is a combination of arclength continuation and branch switching. Briefly speaking, given an initial guess ðu 0 ; λ 0 Þon a branch, arc-length continuation will trace out the remaining part of this branch along the variations of the parameter λ. On the other hand, branch switching will detect bifurcation points along the branch and construct initial solutions on branches emanating from it. Once one solution on each emanating branch is computed, arc-length continuation is applied to complete the remaining part of the new branch.
This combination of arc-length continuation and branch switching is very powerful for computing connected bifurcation diagrams. However, in the presence of geometric imperfections that disconnect branches, it fails to compute these other branches. One approach is to restore the broken symmetry group, find all branches of the now continuous diagram via branch switching, re-introduce the asymmetry, and continue the solutions from the symmetric to the asymmetric state. If the asymmetry is introduced via a parameter in the equations (e.g. asymmetry in the loading), this is straightforward, but if the asymmetry is introduced in other ways (such as in the geometry) this procedure can be very difficult to apply. Farrell et al. [54] remedy this issue by introducing the deflated continuation algorithm to discover disconnected branches from known ones without requiring any detection of the bifurcation point. The algorithm is used in this work for computing the bifurcation diagrams of different stiffener types.
At the heart of the deflated continuation method is the deflation technique [53]. Historically, the deflation technique was first applied to finding distinct solutions to scalar polynomials [79]. Brown and Gearhart [80] then extended this deflation approach to solving systems of nonlinear algebraic equations via the construction of deflation matrices. A more recent study of Farrell et al. [54] extended the deflation technique to the case of infinite-dimensional Banach spaces, appropriate for partial differential equations. In the following, we recall the idea of deflation.
For a fixed parameter λ � , the parameter-dependent problem (4.1) becomes Suppose that (4.2) permits multiple solutions and the Newton iteration converges to a known solution u � . The goal of deflation is to find as many solutions to FðuÞ ¼ 0 in a way that the Newton iteration will never converge to known solutions even with the same initial guess. To this end, a new problem GðuÞ � Mðu; u � ÞFðuÞ ¼ 0 is constructed, where Mðu; u � Þ is a deflation operator and G satisfies the following two properties: 1 GðuÞ ¼ 0 has the same solutions as those of FðuÞ ¼ 0; that is to say, for all u 6 ¼ u � , FðuÞ ¼ 0 ⇔ GðuÞ ¼ 0. 2 For a known solution u � to FðuÞ ¼ 0, G will not converge to u � again; i.e., given any sequence u i →u � , liminf ui→u � jjGðu i Þjj > 0.

The form of Mðu; u � Þ used in this work is the shifted deflation operator
where the pole strength p governs the rate at which the function approaches the introduced singularity, and the shift parameter α ensures that the deflated problem recovers the behaviour of the original problem far from previously found solutions as ku u � k→∞. In our algorithm, the values p ¼ 2 and α ¼ 1 are adopted.
We now give a brief description of the deflated continuation algorithm. The algorithm proceeds by continuation over a range of values of ε. Consider the step in the algorithm going from ε ¼ ε to ε ¼ ε þ . Suppose that n solutions u 1 ; u 2 ; …; u n are known at ε ¼ ε . The step proceeds in two phases. First, each solution u i is continued from ε to ε þ yield u þ i (using arclength, tangent or standard continuation). 3 As each solution u þ i is computed, it is deflated away from the nonlinear problem at ε ¼ ε þ . Once all known solutions have been continued, the search phase of the algorithm begins. Each previous solution u i is used again as initial guess for the nonlinear problem at ε ¼ ε þ ; the deflation operator ensures that the solve will not converge to any of the known solutions u þ i , and hence if Newton's method converges it must converge to a new, unknown solution. Importantly, this unknown solution may lie on a disconnected branch. If an initial guess yields a new branch, the new solution is deflated and the initial guess used repeatedly until failure.
Once all initial guesses from ε have been exhausted, the step completes and the algorithm proceeds to the next step. This is repeated until the continuation parameter reaches a desired target value. The search is applied at all steps, i.e. no a priori knowledge of the location of the disconnected bifurcations is assumed. For more details, including application to standard benchmark cases, see Ref. [54].

Bifurcation analysis of buckling behaviours
In this subsection, the deflated continuation algorithm is applied to investigate the buckling behaviour of three different aircraft stiffeners, i. e., the L-shaped asymmetric profile, the L-shaped symmetric profile and the Z-shaped asymmetric profile (see Fig. 1). We perform a uniaxial compression test along the z-axis, with the compressive force applied to the bottom face of each stiffener. SI units are adopted for all physical quantities in the subsequent experiments.
Throughout the simulations, the boundary data ε is continued with Fig. 2. Displacements (in mm) along the y-axis showing the bifurcation path switch due to a small perturbation load on the L-shaped asymmetric stiffener with a bulb.

Fig. 3. Displacements (in mm) along the x-axis
showing the bifurcation path switch due to a small perturbation load on the Z-shaped asymmetric stiffener. 3 In the problems considered here standard zero-order continuation is sufficient, so we use this. the continuation step Δε ¼ 10 5 . All numerical experiments are based on a continuous piecewise linear discretisation of the displacement function space. For the linearisation, we employ Newton's method with the L 2 line search algorithm of PETSc [64] with relative and absolute tolerances 10 7 . The solve is terminated with failure if convergence is not achieved in 50 nonlinear iterations. At each Newton iteration, the linearised system is solved by GMRES with a V-cycle multigrid preconditioner, where the coarse grid problem is solved by Cholesky factorisation and the additive block Successive Over-Relaxation (SOR) algorithm is used as relaxation [64]. The computations are performed on eight cores, parallelised using MPI.

L-shaped asymmetric profile
In this experiment, the boundary data ε is varied in the range ½0; 0:002�. Fig. 4 illustrates the bifurcation diagram of the functional u 1 ð0:019; 0; 0:025Þ with respect to the parameter ε. For ε≲0:00078 , there is only one solution to problem (2.10); two more solutions appear until ε � 0:00086 , after which there exist at least five solutions. Then from ε � 0:00139, two more branches are found, leading to a total number of seven solutions discovered.
The resulting seven solutions at ε ¼ 0:002 are given in Fig. 5. For a better connection with the bifurcation diagram in Fig. 4, we point out that the first and second deformed profiles in Fig. 5 correspond to the lowest and uppermost branch.
Furthermore, we compute the stability through calculating the inertia of the Hessian matrix of the energy function ΦðuÞ with a Cholesky factorisation [81], Section 16.2]. This reveals that the first two buckling profiles in Fig. 5 are stable (i.e., the Hessian matrix is positive definite) while the remaining five modes (with a nonzero number of negative eigenvalues, making the Hessian matrix indefinite) are unstable. These five unstable buckling profiles can be easily perturbed.
From Fig. 4, it is noticeable that there exist disconnected branches even though the body force and the traction are zero in the model. This is due to the non-symmetric geometry of the aircraft stiffener, making it easier to buckle outwards than inwards.
Additionally, we plot the average stress over the bottom face computed by (2.6) in Fig. 4. It is shown that for sufficiently small deformations, it is proportional to the displacement, as expected from Hooke's law. For relatively large deformations, their relationship becomes nonlinear. Notice that the yielding stress of Aluminum is 70 MPa [82] and our obtained stress is about the level of 1000 MPa, as can be seen from Fig. 4. Other mathematical models [83] that include plasticity should therefore be considered in the future.

Remark 4.1.
One might wonder about the utility of identifying the unstable buckling modes presented in Fig. 5 above and Figs. 7 and 9 below, as only stable solutions can be physically observed in experiments. However, unstable solutions provide important information about the energy barrier that the system must overcome to switch from one stable solution to another. Of all possible paths in the energy landscape, the one with lowest energy cost will go through one of those unstable solutions (a mountain pass). This intuitive statement is formalised by the so-called Mountain Pass Theorem, see Ref. [84], Section 8.5]. Therefore, knowledge of the unstable modes gives knowledge of the energetic stability of the different local minimisers.

L-shaped symmetric stiffener
We conduct similar numerical experiments for the L-shaped symmetric stiffener which possesses a geometric symmetry due to the absence of the bulb (see Fig. 1). In our preliminary experiments, we observe many more branches in the bifurcation diagram for ε 2 ½0; 0:002�. To make a clear bifurcation diagram, we instead illustrate the case of varying ε in ½0; 0:0007�.
The bifurcation diagram of the functional u 1 ð0:019; 0; 0:025Þ is shown in Fig. 6. We first observe that when ε≲0:00031, there exists only one solution. The system then undergoes a pitchfork bifurcation with three solutions until ε � 0:00047, after which it presents five solutions. Around the point of ε � 0:00065, it starts to buckle in seven different modes with two new disconnected branches.
One expects a connected bifurcation diagram for this stiffener profile because of its geometric symmetry. However, Fig. 6 reveals a disconnected bifurcation around ε � 0:00065. This is the correct diagram for this discrete problem, and the reason is subtle: while the mesh is almost perfectly symmetric, there is a slight asymmetry around the centre web. In general, exactly preserving the continuous symmetries of the geometry during mesh generation is very difficult. Even small perturbations to the symmetry can lead to a (discrete) disconnected bifurcation diagram that may not be easily captured by conventional arc-length continuation and branch switching algorithms. The deflated continuation algorithm we used helps us capture the relevant branches without needing to enforce symmetry of the mesh, improving the flexibility of the computations.
In this numerical experiment, we have found seven buckling modes in total and they are all illustrated in Fig. 7 in pairing order. There are three Z 2 -symmetric pairs of modes, as well as the single Z 2 -symmetric compressed state. Additionally, the stability of each buckling profile is indicated in Fig. 7. We can see that only the first two profiles are stable while the remaining five buckling profiles can be easily perturbed.
The average stress over the bottom face computed by (2.6) is plotted in Fig. 6. We can observe a linear stress-strain relation for small ε corresponding to Hooke's law and then this relation becomes nonlinear for larger deformations. As before, a more physically realistic model should incorporate plasticity.

Z-shaped asymmetric stiffener
For the Z-shaped asymmetric stiffener profile, its bifurcation diagram is shown in Fig. 8. To keep the number of solutions considered manageable, we consider ε 2 ½0; 0:0027�. The disconnection of the bifurcation comes again from the asymmetry of the domain, similar to the case of the L-shaped asymmetric stiffener. It can be seen that the diagram starts to bifurcate at ε � 0:00191, obtaining three solutions, and approximately at ε ¼ 0:00203, five branches appear until ε � 0:00249 where four more solutions are found. Consequently, there are nine branches in total and Fig. 9 illustrates these buckling profiles at ε ¼ 0:0027. Essentially, four pairs of buckling modes have been discovered, along with the neutrally compressed state.
We also point out that the uppermost and lowest branch in Fig. 8 correspond to the first and the second buckling profiles in Fig. 9. This implies that the Z-shaped asymmetric stiffener is easier to buckle upwards rather than downwards.
Regarding the stability of each buckling profile, it is shown that only the first two buckling profiles in Fig. 9 are stable (i.e., the Hessian matrix is positive definite). The remaining seven unstable solutions in Fig. 9 can be easily perturbed.
The average stress over the bottom face, computed by (2.6), is plotted in Fig. 8. The linear stress-strain relation for small ε is also observed, which again verifies Hooke's law, and again indicates that plasticity should be considered to achieve more physically realistic results.

Conclusions and future work
The aim of this paper has been to investigate multiple post-buckling bifurcation paths for three types of thin-walled stiffeners. Deflated continuation allows for the effective capturing of more bifurcation paths, especially disconnected ones. For each bifurcation path, its stability is calculated, allowing more reliable analysis of the buckling profiles. In future work, deflation should be applied to larger structures, including shell models and models incorporating plasticity. This development will enable the investigation of multiple bifurcation paths expected in highly optimised large structures such as aircraft wings, in which the skin pockets in different regions usually buckle at similar load levels. In these cases, the computation of the disconnected bifurcation diagram can generate a robust basis for comparison with experimental results, in the sense that the experimental results should correspond to one of the bifurcation branches. The authors also suggest applying the deflation technique to the design of stiffened panels which are prone to mode jumps, aiming to achieve designs that are free from this undesirable post-buckling behaviour.  We denote the left-hand-side bilinear form by Aðu r ; v; δuÞ and the L 2 -norm over the Dirichlet boundary Γ D by k ⋅k 0;ΓD . To ensure the solvability of the above Newton system, the major issue is to prove that with the current approximation u r 2 V, Aðu r ; w; wÞ > 0 for all nonzero w 2 V. In the following, we will ignore the superscript r for notational simplicity.
Note that and Young's inequality (ab � εa 2 þ b 2 4ε ). One should notice that the constant C 2 in the inverse inequality (A.2) scales like the bulk modulus K denoted by Since the parameters E ¼ 69 GPa, ν ¼ 0:334 are chosen for aluminium stiffeners [72,73], the bulk modulus is K � 70 GPa by (A.3). In addition, the material we considered in this work is compressible [66] as ν < 0:5.