One-variable Optimization of the CRM-like Box Beam¶

In our last notebook we developed a straight box beam approximation of the CRM wingbox, which we called CRM-like box beam, and we studied its nonlinear structural response under a concentrated tip load. In this notebook we are going to use the OpenMDAO framework to attempt a proof-of-conpcept optimization of such box beam model including a nonlinear structural stability constraint. In this proof-of-concept optimization we are going to minimize the mass of the structure using a single design variable, that is to say the wall thickness, which will be the same for all structural elements. After recalling the CRM-like box beam geometry from our last notebook, we will define the optimization problem in a way that the initial structure is optimal in linear buckling terms, to allow an easy comparison. Then we'll implement the optimization problem with the OpenMDAO library, set up the numerical model and carry the optimization taking a look at how some nonlinear analysis parameters affect the results.

• Definition of the CRM-like box beam model
• Formulation of the optimization problem
• OpenMDAO implementation
• Setup of the numerical model
• Optimization
• Conclusion

Definition of the CRM-like box beam model ¶

---

We recall the geometrical parameters that we derived in our last notebook.

We also recall the target element length that we obtained from the mesh convergence study.

Then we recall the material properties. This time we also consider the yield strength to enforce material strength constraint during the optimization. All material properties are taken from the paper by Brooks et al. (2018).

Formulation of the optimization problem ¶

---

In this notebook we want to attempt a proof-of-concept optimization with a nonlinear structural stability constraint. Given the proof-of-concept nature of our attempt, we consider only one design variable, the wall thickness $t$, which we consider to be the same for all the structural elements. We define our optimization objective as the minimization of structural mass, and we bound the wall thickness between a minimum and a maximum value.

What load magnitude should we apply to our box beam? The idea is that we should be able to compare the structure optimized with the nonlinear structural stability constraint with a baseline design that is optimal in linear buckling terms.

With the wall thickness being our only design variable, we can immediately notice that there is a bijective relationship between the wall thickness and the linear buckling load of the box beam, which means that for each value of the wall thickness there is a corresponding value of the linear buckling load. If we choose the design load to be equal to the linear buckling load of the initial structure with wall thickness $t_0$, we can state that the initial structure is already optimal in linear buckling terms. In fact, if we were to decrease the wall thickness, the linear buckling load would become lower than the design load, making the structure unfeasible.

As a consequence, by setting up our optimization problem with only one design variable and by applying a design load equal to the linear buckling load of the initial structure, we can assess the efficacy of our nonlinear structural stability constraint by simply comparing the optimized structure to the initial design.

Let's recall the linear buckling load of the CRM-like box beam that we found in our last notebook.

Now we need to define the constraints of our optimization. We are going to consider the following constraints:

• the structure must be stable at every equilibrium point along the equilibrium path;
• all deformations must be elastic, i.e. the material must not yield.

The first constraint corresponds to our nonlinear structural stability constraint, and we are going to enforce it by constraining the $N_\lambda$ smallest-magnitude eigenvalues of the tangent stiffness matrix to be positive for every converged increment of the analysis. The second constraint corresponds to a more traditional material strength constraint, and we are going to enforce it by constraining the von Mises stress of all elements to remain below the yield strength.

Our optimization problem can be stated as it follows: $$ \begin{aligned} \text{minimize} & & & m\\ \text{by varying} & & & t_\mathrm{min} \leq t \leq t_\mathrm{max}\\ \text{subject to} & & & g_{\lambda,\,ij}=-\lambda_{ij}<0 & & i=1, \dots, N_i,\; j=1, \dots, N_\lambda\\ & & & g_{\sigma,\,k}=\sigma_k-\sigma_\mathrm{max}<0 & & k=1, \dots, N_e, \end{aligned} $$

where $t_{min}$ and $t_{max}$ are the thickness bounds, $\lambda_{ij}$ is the $j$-th tangent stiffness matrix eigenvalue of the $i$-th increment, $\sigma_k$ is the von Mises stress of the $k$-th element, and $\sigma_{max}$ is the yield strength of the material.

Instead of imposing the constraints on each eigenvalue and on each element, we are going to aggregate them using Kreisselmeier–Steinhauser (KS) functions (Martins & Ning, 2021), so that we have one single value for each constraint function:

$$KS_\lambda=\max_{i,j}\left(g_{\lambda,\,ij}\right)+\frac{1}{\rho}\ln\left(\sum_{i=1}^{N_i}\sum_{j=1}^{N_\lambda}\exp\left(\rho\left(g_{\lambda,\,ij}-\max_{i,j}\left(g_{\lambda,\,ij}\right)\right)\right)\right),$$$$KS_\sigma=\max_k\left(g_{\sigma,\,k}\right)+\frac{1}{\rho}\ln\left(\sum_{k=1}^{N_e}\exp\left(\rho\left(g_{\sigma,\,k}-\max_k\left(g_{\sigma,\,k}\right)\right)\right)\right).$$

Furthermore, we need to include an additional constraint for the applied load at the end of the nonlinear analysis. In fact, even if the nonlinear analysis fails to converge to the prescribed applied load, Nastran will still produce valid output files that will be fed to the optimizer. When the applied load at the end of the nonlinear analysis $P_\text{end}$ is lower than the prescribed design load $P_\text{design}$, the output file will provide misleading values for the constraint functions. As a consequence, the following constraint is defined:

$$g_P=\frac{P_\text{end}}{P_\text{design}}\geq 0.99,$$

which means that the applied load at the end of the nonlinear analysis must be at least 99% of the prescribed design load.

OpenMDAO implementation ¶

---

To implement our optimization, we are going to use the OpenMDAO libary, which is an open-source high-performance computing architecture for multidisciplinary analsyis and optimization. The smallest unit of computational work within OpenMDAO is the so-called Component, where the actual numerical calculations happen. Each component has input and output variables that can be connected to chain together multiple computations. In our case, we only have one kind of computation happening, that is to say the nonlinear analysis under the design load. This computation has to take as input the wall thickness of our box beam, which is our only design variable, and calculate our objective and constraints.

Let's go ahead and define this component, that we are going to call NastranSolver. We define the component such that it takes the following variables as input options:

• the BDF object, representing the Nastran model;
• the input name of the Nastran run;
• the flag to run the Nastran analysis.

These options represent "static inputs" of the component, defining quantities or objects that do not change during the optimization. The real input of the component, or in other words the quantity that does change during the optimization, is the wall thickness of the box beam, and we define it accordingly.

Finally, we define the following component outputs:

• the mass of the box beam;
• the value of the nonlinear structural stability constraint;
• the value of the material strength constraint;
• the value of the applied load constraint.

Then we define the compute method of the component such that it changes the wall thickness of the box beam in the BDF object, calculates the mass of the box beam, runs the nonlinear analysis, and calculates the constraints values.

Every OpenMDAO model is constructed using instances of the Group superclass, forming a hierarchy that arranges the components. In our case the hierarchy is very simple, as we only need a single group holding our NastranSolver. We then define the NastranGroup, which takes the same options of the NastranSolver component.

With the NastranSolver and the NastranGroup classes we have defined the computational core of our optimization problem. We can now move on to the setup of the numerical model.

Setup of the numerical model ¶

---

To define an instance of NastranGroup and set up the actual optimization problem, we first need to define the input BDF object. We do this by generating a mesh with the function mesh_box_beam_reinforced_with_ribs_and_stiffeners and then passing it to the function create_base_bdf_input, both belonging to the box_beam_utils module.

This gives us a base BDF object including material properties, nodes, elements, boundary conditions and output files defaults. We proceed with adding the command to request the output of the nonlinear stresses.

Analogously to our last notebook, we apply a vertical force at the tip of the CRM-like box beam, setting the magnitude of the force equal to the value of the design load. Successively, we create the corresponding static subcase.

Now we need to setup the nonlinear analysis with the arc-length method. What parameters should we use? The key requirements are to trace the equilibrium path of the structure in a robust way, or in other words avoiding to jump away from it, and to do so in a computationally efficient manner. As a consequence, we cannot choose to employ very fine arc-length increments from the start to the end of the analysis, as this would result in exessively long computational times.

In our last notebook, we observed that by using load and displacement error functions, $\mathrm{MAXITER}=5$ and remaining default parameters, all the equilibrium points found by the nonlinear analysis lied on the natural equilibrium path. This analysis took only about 200 seconds, so only a fraction of the 3.2 hours taken by the analysis with very fine arc-length increments. As a consequence, we might be tempted to use the mentioned parameters for the nonlinear analyses of our optimization problem.

However, the analysis with load and displacement error functions, $\mathrm{MAXITER}=5$ and remaining default parameters was not able to capture the unstable part of the equilibrium path and thus the snap-through behavior of the structure. Capturing the unstable part of the equilibrium path is an essential requirement for the nonlinear analyses of our optimization problem, in order to enforce the nonlinear structural stability constraint. As a consequence, we need to refine this set of parameters with the aim of capturing the snap-through behavior of the structure.

To start, we read the op2 file of the analysis employing the mentioned parameters and of that employing very fine arc-length increments. We then plot the 3D equilibrium diagram showing the $z$-displacement at the root of the CRM-like box beam, the $z$-displacement at the tip and the applied load. We do this by using the function plot_3d_load_displacements_stability from the plot_utils module to highlight the stability along the equilibrium paths.

As described earlier, we can see that the analysis with load and displacement error functions and $\mathrm{MAXITER}=5$ is not able to capture the unstable part of the equilibrium path.

All the remaining parameters for this analysis are set to their default values, including the desired number of iterations, $I_d$, which is set to $12$. This parameter controls the ratio of the step size between successive iterations:

$$\frac{\Delta l_{new}}{\Delta l_{old}}=\sqrt{\frac{I_d}{I_{max}}},$$

where $I_{max}$ is the number of iterations required to achieve convergence for the last step.

We notice that by choosing $\mathrm{MAXITER}=5$ and $I_d=12$ we are effectively imposing a step size ratio equal or larger than $\sqrt{2.4}$ ($\sqrt{12/5}$). This means that the arc-length solver will not be able to take steps smaller than $\sqrt{2.4}$ times the step size of the last converged step, unless using increment bisection. It then makes sense to set $I_d$ to a value equal to $5$, so that the solver can use step size ratios between $1$ and $\sqrt{5}$, corresponding respectively to $I_{max}=5$ and $1$.

We set up the nonlinear analysis with load and displacement error functions, $\mathrm{MAXITER}=I_d=5$, apply twice the design load, run the analysis, and visualize the results.

The analysis with $I_d=5$ fails to converge to the prescribed applied load because it reaches the default maximum number of arc-length increments, $20$. At the same time, we observe a jump from the natural path to a complementary path at the fourth increment. We inspect the convergence results of the analysis at this increment in the f06 file.

                        - - - CONVERGENCE  FACTORS - - -                      - - - LINE SEARCH DATA - - -
0ITERATION      EUI         EPI         EWI      LAMBDA       DLMAG       FACTOR     E-FIRST     E-FINAL   NQNV  NLS  ENIC  NDV  MDV

         3  3.2328E-03  2.9659E-03  3.7872E-05  1.0000E-01  4.8705E-02  1.0000E+00  7.4023E-02  7.4023E-02    0    0          0    1
 *** USER INFORMATION MESSAGE 6186 (NCONVG)
     *** SOLUTION HAS CONVERGED ***
         SUBID    1     LOOPID    4     LOAD STEP    0.586     LOAD FACTOR 0.58599681

We can see that the increment converges at the third iteration, with a displacement error of $3.2\cdot 10^{-3}$ and a load error of $3.0\cdot 10^{-3}$.

We then lower both convergence tolerances to $10^{-3}$ and we increase the maximum number of increments to $100$. We also lower the minimum allowable arc-length adjustment ratio to $10^{-2}$, as the bisections of the increments in the last analysis resulted in arc-length adjustment ratios in the order of $10^{-2}$.

In this case, the analysis fails to converge to the prescribed applied load because it reaches the maximum number of bisections. Also here, we observe a jump from the natural path to a complementary path at the fourth increment. Once again we inspect the convergence results of the analysis at this increment in the f06 file.

                        - - - CONVERGENCE  FACTORS - - -                      - - - LINE SEARCH DATA - - -
0ITERATION      EUI         EPI         EWI      LAMBDA       DLMAG       FACTOR     E-FIRST     E-FINAL   NQNV  NLS  ENIC  NDV  MDV

         4  1.6871E-04  3.1931E-05  7.3704E-09  1.0000E-01  2.0420E-04  1.0000E+00  0.0000E+00  0.0000E+00    0    0          0    1
 *** USER INFORMATION MESSAGE 6186 (NCONVG)
     *** SOLUTION HAS CONVERGED ***
         SUBID    1     LOOPID    4     LOAD STEP    0.611     LOAD FACTOR 0.61064795

The increment converges at the fourth iteration, with a displacement error of $1.7\cdot 10^{-4}$ and a load error of $3.2\cdot 10^{-5}$.

We modify again the nonlinear analysis parameters, setting the displacement tolerance to $10^{-4}$, the load tolerance to $10^{-5}$, the maximum number of bisections to $20$ and the maximum number of iterations and the desired number of iterations to $4$. We run the analysis and visualize the results.

Eureka! We've succeded in converging to the perscribed applied load and in capturing the unstable part of the natural equilibrium path. We can observe that the selected nonlinear analyses parameters have allowed the arc-length solver to take relatively large increments at the beginning of the analysis and to refine them as the structure approaches the load limit point and the unstable part of the equilibrium path. In this way we have achieved a good compromise between computational efficiency and robustness, employing approximately 25 minutes to run the nonlinear analysis.

Can we find a way to further reduce the computational time? A possible solution is to consider an alternative approach to our current full Newton iteration method, where the tangent stiffness matrix is updated at every iteration. A popular alternative is the so-called modified Newton iteration method, where the tangent stiffness matrix is updated only at the beginning of the increment and then kept constant. This approach is less computationally expensive than the full Newton iteration method, but it is also less robust, because the solver does not use the true Jacobian of the system. The implication is that convergence can be slowed down and in extreme cases it can be prevented at all. Consequently, we need to test the modified Newton iteration method with the current nonlinear analysis parameters to verify whether it can provide a computational time reduction without compromising the robustness of the analysis.

The closest implementation of the modified Newton iteration method in Nastran can be obtained using $\mathrm{KSTEP}=\mathrm{AUTO}$. With this option, the solver automatically selects whether to update the stiffness matrix at each iteration or use the old matrix based on convergence rates. In practice, this means that the solver always tries not to update the stiffness matrix, unless it estimates that the number of iterations to converge exceeds $\mathrm{MAXITER}$ or that the time required for convergence with current stiffness exceeds the estimated time required for convergence with updated stiffness. The tangent stiffness matrix is always updated when the increment converges, so that it can be used for the next increment.

We set $\mathrm{KSTEP}$ to $\mathrm{AUTO}$ and we also increase the maximum number of increments to 200 and decrease the minimum allowable arc-length adjustment ratio to $10^{-5}$, as the bisections of the increments in the last analysis resulted in arc-length adjustment ratios in the order of $10^{-5}$.

The analysis fails to converge to the prescribed load because it reached the maximum number of bisections. Looking at the load-displacement diagram, we can see that the analysis stops a few increments after reaching the load limit point. Inspecting the f06 file, after the last converged increment we find several instances of no root case being encountered:

                        - - - CONVERGENCE  FACTORS - - -                      - - - LINE SEARCH DATA - - -
0ITERATION      EUI         EPI         EWI      LAMBDA       DLMAG       FACTOR     E-FIRST     E-FINAL   NQNV  NLS  ENIC  NDV  MDV

         1  2.3212E-04  3.2444E-05  3.1665E-10  1.0000E-01  8.1316E-04  1.0000E+00  0.0000E+00  0.0000E+00    0    0          0    1
 *** USER INFORMATION MESSAGE 6192 (NLROOT)
     *** NO ROOT CASE HAS BEEN ENCOUNTERED ***
 *** USER INFORMATION MESSAGE 6187 (NLITER)
     *** BISECTION METHOD IS NOW ACTIVATED ***
 *** USER INFORMATION MESSAGE 6188 (NLITER)
     ***  ARC LENGTH RATIO IS 6.103516D-06
         1  1.1606E-04  2.0971E-05  1.2499E-10  1.0000E-01  7.1734E-04  1.0000E+00  0.0000E+00  0.0000E+00    0    0          0    1
 *** USER INFORMATION MESSAGE 6192 (NLROOT)
     *** NO ROOT CASE HAS BEEN ENCOUNTERED ***
 *** USER INFORMATION MESSAGE 6187 (NLITER)
     *** BISECTION METHOD IS NOW ACTIVATED ***
 *** USER INFORMATION MESSAGE 6188 (NLITER)
     ***  ARC LENGTH RATIO IS 3.051758D-06
         1  5.8029E-05  1.5320E-05  5.4208E-11  1.0000E-01  6.6947E-04  1.0000E+00  0.0000E+00  0.0000E+00    0    0          0    1
 *** USER INFORMATION MESSAGE 6192 (NLROOT)
     *** NO ROOT CASE HAS BEEN ENCOUNTERED ***
 *** USER INFORMATION MESSAGE 6187 (NLITER)
     *** BISECTION METHOD IS NOW ACTIVATED ***
 *** USER INFORMATION MESSAGE 6188 (NLITER)
     ***  ARC LENGTH RATIO IS 1.525879D-06
         1  2.9014E-05  1.2531E-05  2.5045E-11  1.0000E-01  6.4554E-04  1.0000E+00  0.0000E+00  0.0000E+00    0    0          0    1
 *** USER INFORMATION MESSAGE 6192 (NLROOT)
     *** NO ROOT CASE HAS BEEN ENCOUNTERED ***
 *** USER INFORMATION MESSAGE 6187 (NLITER)
     *** BISECTION METHOD IS NOW ACTIVATED ***
 *** USER INFORMATION MESSAGE 6188 (NLITER)
     ***  ARC LENGTH RATIO IS 7.629395D-07
         1  1.4507E-05  1.1155E-05  1.2011E-11  1.0000E-01  6.3358E-04  1.0000E+00  0.0000E+00  0.0000E+00    0    0          0    1
 *** USER INFORMATION MESSAGE 6192 (NLROOT)
     *** NO ROOT CASE HAS BEEN ENCOUNTERED ***
 *** USER INFORMATION MESSAGE 6193 (NLITER)
     *** MAXIMUM NUMBER OF BISECTIONS OR MINIMUM ARC LENGTH HAS BEEN REACHED.
 ^^^     
 *** USER INFORMATION MESSAGE 9052 (SUBDMAP NLSTATIC)
     NONLINEAR STATIC ANALYSIS COMPLETED.

This results suggest that the arc-length solver is struggling to find solutions on the unstable part of the equilibrium path without updating the stiffness matrix.

To solve this issue, we can try to set $\mathrm{KMETHOD}$ to $\mathrm{SEMI}$. With this option, the solver carries out the first iteration of the new increment, updates the stiffness matrix, and resumes with $\mathrm{KMETHOD}=\mathrm{AUTO}$. By always updating the stiffness matrix after the first increment, we expect the robustness of the solver to increase. Let's test this on our case.

Great, the analysis has successfully converged to the prescribed load! However, looking at the computational time, it seems that we did not achieve in our goal of reducing it. In fact our last analysis took about 80 seconds longer than the analysis with the full Newton method. Consequently, we choose to use the latter to run the optimization.

To assess the accuracy of the analysis with the efficient set of parameters, we calculate the percentage difference of the applied load at the first limit point with respect to the reference analysis with very fine arc-length increments.

The negligible percentage difference confirms that by choosing the more efficient set of parameters accuracy is not lost with respect to the reference analysis.

Before moving on to the optimization, is there anything else that we can do to reduce the computational time of the nonlinear analyses? Yes there is! We can use the parallel computing capabilities of SOL 106. Since it is quite an old solution sequence, SOL 106 does not have great parallel computing capabilities, but it can still parallelize matrix decomposition and multiplication. Better than nothing!

To use the parallel computing capabilities of SOL 106 we need to specify NASTRAN PARALLEL=N in the system command section of the input file, where N is the number of cores. How many cores should we use? Well, of course it depends on the machine that we have! On my machine I have 12 logical processors, meaning that we can use up to 12 cores. To test what number of cores gives the smallest computational time, we run the analysis with 2, 4, 8 and 12 cores, and look at the final computational times.

The analysis with 4 cores is the one that results in the lowest computational time, so we add NASTRAN PARALLEL=4 to the system command section.

Optimization ¶

---

We are finally ready to set up and run our first optimization. We start by restoring the mangnitude of the applied load to that of the design load and by reducing the maximum number of increments to 100, as the analysis employing the efficient set of parameters employed 63 iterations to converge.

Then we set up the nonlinear analysis to calculate the 20 smallest eigenvalues of the tangent stiffness matrix for every converged increment, by calling the function set_up_sol_106_with_kllrh_eigenvalues. This function does the same set up work of run_sol_106_with_kllrh_eigenvalues but it does not run the analysis. To further reduce the computational time of the optimization, this time we add the option dmap_option='stop'. This option includes in the bdf file a different version of the DMAP for the calculation of the eigenvalues of the tangent stiffness matrix, which stops the analysis when the first negative eigenvalue is found. In this way we can avoid running the analysis for increments that are beyond the loss of stability, since these have no interest for the optimization.

To run our optimization with OpenMDAO we need to define an instance of the Problem class. This class constitutes the top-level container for the model hierarchy, and it is used to set up the optimization problem. Let's create an instance of this class using the previously defined NastranGroup as the model argument. We then call the setup method to get the data structures of the framework in place.

A good practice in optimization is to normalize objective, design variables and constraints with respect to some reference values. A very popular choice for this normalization is to use the values corresponding to the initial design. For this reason, we assign the initial value of wall thickness to the design variable, run our Problem object with the method run_model, and record the values of the objective and constraint functions.

Now we set up the optimization problem in the follwowing steps:

• define the default input values of the NastranSolver component;
• select COBYLA as the optimization algorithm;
• define the design variable, objective and constraints functions together with their bounds and scaling factors.

For the wall thickness we select an bounding interval between 1 and 20 mm.

Before running the optimization, we are going to create a recorder and attach it to our Problem object. This will allow us to save the values of the design variable, objective and constraints at each iteration of the optimization.

Finally, we can run the optimization calling the method run_driver of our Problem object.

Let's plot the convergence history of the optimization. To do this, we define the function plot_optimization_history, that takes as input the path to the recorder file, plots the convergence history of the design variable, constraints and objective, and returns a dictionary with their histories.

The optimization history shows that the last design point is actually unfeasible, as $KS_\lambda$ has a positive value. We consider the last feasible design point and find the final value of wall thickness and of the constraints. We also calculate the percentage mass reduction with respect to the initial structure.

From the above results, it seems that the optimizer managed to reduce the wall thickness only slighlty, obtaining a meager 1% mass reduction. The margin of stability is also reduced only slightly with respect to that of the initial structure, suggesting that the structure has still much margin for mass reduction. What is happening?

Let's plot the eigenvalue and the load-displacement diagram of the structure at the last point of the optimization, the unfeasible one, to inspect the response obtained from the nonlinear analysis.

The load-displacement diagram shows that the final unstable equilibrium point is characterized by a larger root displacement magnitude compared to the other equilibrium points, but with a negative value. This is in contrast with the results of the initial structure, where there is indeed an increase in the root displacement corresponding to the broken supercritical pitchfork, but the displacement remains positive. This suggests that the arc-length solver might have accidentally jumped on a complementary path. We then inspect the convergence of the last increment in the f06 file.

                        - - - CONVERGENCE  FACTORS - - -                      - - - LINE SEARCH DATA - - -
0ITERATION      EUI         EPI         EWI      LAMBDA       DLMAG       FACTOR     E-FIRST     E-FINAL   NQNV  NLS  ENIC  NDV  MDV

         4  9.9904E-05  3.1832E-06  1.1991E-09  1.0000E-01  4.8751E-05  1.0000E+00  0.0000E+00  0.0000E+00    0    0          0    1
 *** USER INFORMATION MESSAGE 6186 (NCONVG)
     *** SOLUTION HAS CONVERGED ***
         SUBID    1     LOOPID    7     LOAD STEP    1.000     LOAD FACTOR 1.00000000

The increment converges at the fourth iteration, with a displacement error of $9.9\cdot 10^{-5}$ and a load error of $3.2\cdot 10^{-6}$. To verify whether the solver has jumped to a complementary path we can repeat the analysis using tighter set of parameters to enforce convergence. Looking at the convergence results of the last increment, we can do this in two way:

• reduce $\mathrm{MAXITER}$ from 4 to 3 keeping the same convergence tolerances;
• or reduce the convergence tolerances to $\mathrm{EPSU}=10^{-5}$ and $\mathrm{EPSP}=10^{-6}$ keeping $\mathrm{MAXITER}=4$.

Let's try with the second option, with the aim of keeping $\mathrm{MAXITER}=4$ and providing the solver with more room to increase the increment where the equilibrium path is not extremely nonlinear. We run the analysis and we plot the eigenvalue and the load-displacement diagrams.

The above results suggest indeed that in the previous analysis the solver jumped to a complementary path because the convergence tolerances were too loose for the analyzed case. In fact, in the last analysis, the solver reduced the arc-length increments drastically as the applied load approached to $P/P_\text{desgin}=1$.

Let's repeat the optimization with these values of convergence tolerances.

Once again, the optimization ends on an unfeasible point, so we look at the design variable, constraints and objective of the last feasible point.

The results show an analogous situation to what we found after the first optimization. We then plot the eigenvalue and load-displacement diagrams to check whether the arc-length solver has jumped to a complementary path in the last analysis.

The above results are totally analogous to what we found earlier, confirming that also in this case, despite the use of tighter convergence tolerances, the solver jumped to a complementary path because the convergence tolerances were still too loose.

Consequently, we have to consider the other option mentioned earlier, that is to say lowering $\mathrm{MAXITER}$ from 4 to 3. However, what value of desired number of iterations should we use? If we lower also $I_d$ to 3, then the arc-length solver will have very small room to increase the arc-length increments, as the only case when this will happen is when the increment converges in 2 iterations. This can cause the arc-length increments to remain small during the analsis, entailing an increase in the computational time. At the same time, if we keep $I_d=4$ to give more room to the solver to take larger increments, the increments will always be increased. This can entail an excessive increase of the increments and cause the solver to take more bisections.

To make an inoformed choice between $I_d=3$ and $4$, we run two different analyses on the initial structure. We apply twice the design load and stop the analysis at the first negative eigenvalue, to deal with a representative case of the optimization. For each analysis we look at the computational time and we plot the eigenvalue and the load-displacement diagrams.

The analysis with $I_d=4$ takes about 300 seconds less than that with $I_d=3$. In the load-displacement diagrams we can observe that this reduction in computational time is due to the solver taking less increments for the analysis with $I_d=4$. Inspecting the f06 files, we find that, to our surprise, the analysis with $I_d=4$ takes also less bisections. It is unclear why this happens.

Given the above results, we choose to set $I_d=4$ for the optimization. We then run the optimization with $\mathrm{MAXITER}=3$, $I_d=4$, and the old convergence tolerances.

Also in this case the optimization ends with an unfeasible design, so we look for the last feasible point and print the values of the design variable, contraints and objective.

This time we find a significantly larger reduction of the mass and of the margin of stability, suggesting that optimization did not end because the arc-length solver jumped on a complementary path, rather because it was not possible to find a better design point within the optimizer's tolerances.

To verify that the solver did not jump on a complementary path, we plot the eigenvalue and load-displacement diagrams of the last design, which is unfeasible.

We see that the arc-length solver successfully follows the equilibrium path through the broken supercritical pitchfork without jumping, until it finds the limit point.

Let's run the analysis of the last feasible point.

Before plotting the load-displacement diagram of the optimized structure, let's calculate the linear buckling load to show the glass ceiling of linear buckling.

The above results show that the linear buckling load of the optimized structure is only 71% of the design load, demonstrating the idea of the glass ceiling of linear buckling. In fact, using linear structural stability methods, the optimized sturcture would not be deemed capable of withstanding the design load, whereas the nonlinear analysis suggests that the structure can carry the load without losing stability or failing in terms of material strength. Morover, we observe that the critical buckling mode of the optimized structure is totally analogous to that of the initial structure.

We further demonstrate the idea of the glass ceiling of linear buckling by plotting the load-displacement diagrams of the node at the root and at the tip of the CRM-like box beam, comparing the initial sturcture with the optimized structure.

The above plots show that the optimizer has successfully reduced the mass of the structure by decreasing the margin of structural stability at the design load so that the loss of stability occurs just above such load. The green line shows the load where the optimized structure would be considered to fail for buckling if analyzed with linear methods.

The root displacement diagram shows a large displacement for the optimized structure when compared to the initial one, in an analogous way as it was observed for the initial structure across the broken supercritical pitchfork. As far as the tip displacement diagram is concerned, we can observe that the optimized structure has a globally more compliant response, as it can be expected given the reduced wall thickness.

Given the large root displacement observed for the optimized structure, it is interesting to compare the deformation of the two structures at the design load.

We can clearly observe that the deformation of the initial structure resembles the one of a beam under a linear deformation field, while the optimized structure shows a buckling-like pattern over the root region of the top skin.

Conclusion ¶

---

In this notebook we have demonstrated how to set up and run a proof-of-concept optimization with a nonlinear structural stability constraint.

Initially, we have looked for a computationally efficient set of parameters for the nonlinear analysis, which allows to capture the unstable part of the equilibrium path of the structure without using very small arc-length increments. Successively, we have set up the optimization problem with OpenMDAO, defining the design variable, objective and constraints functions, and we have run the optimization with COBYLA. We further refined the parameters of the nonlinear analysis as we found that the arc-length solver was jumping to complementary paths during the analyses of the optimization, preventing the achievement of significative mass reductions. In this effort, we found out that reducing the maximum number of iterations for each arc-length increment to 3 was the most robust solution.

With the final set of parameters, we have run a final optimization, which resulted in a 10.9% mass reduction with respect to the initial design, which was optimal in linear buckling terms. In other words, the optimized structure resulted in a linear buckling load that is only 71% of the design load, demonstrating the idea of the glass ceiling of linear buckling. The mass reduction was achieved by reducing the margin of structural stability at the design load, so that the loss of stability occurs just above such load. The optimized structure shows a globally more compliant response, as it can be expected given the reduced wall thickness.

In our next notebook we are going to increase the complexity of our optimization problem by considering two design variables defining a linear variation of the wall thickness from root to tip.