In this example we study a more challenging beam problem: The non-axisymmetric buckling of a thin-walled elastic ring, loaded by a spatially-constant external pressure, . For sufficiently small positive (or arbitrarily large negative) values of the ring deforms axisymmetrically and in this mode it is very stiff, implying that large changes in pressure are required to change its radius. However, if exceeds a critical threshold, the axisymmetric configuration becomes unstable, causing the ring to buckle non-axisymmetrically. Once buckled, the ring is much more flexible and small changes in are sufficient to induce dramatic changes in its shape.
The rapid change in stiffness following the buckling makes it difficult to compute strongly buckled solutions by continuation methods as the solution computed for a certain value of may represent a poor approximation of the solution at a slightly larger pressure. In extreme cases, this can cause the Newton method (whose convergence relies on the provision of good initial guesses for the unknowns) to diverge.
We will demonstrate the use of so-called "displacement control" techniques to overcome these problems. Displacement-control techniques are useful in problems in which some a-priori knowledge about the expected displacement field allows us to re-formulate the problem. Rather than prescribing the load on the elastic solid and computing the displacement field, we prescribe the displacement of a carefully-selected material "control" point and regard the load required to achieve this displacement as an unknown.
We wish to compute the deformation of a linearly-elastic, circular ring of undeformed radius and wall thickness , subject to a spatially-constant external pressure . We choose the undeformed radius as the lengthscale for the non-dimensionalisation of the problem and assume that , justifying the use of beam theory. As discussed in the previous example we scale the pressure on the ring's effective Young's modulus, , where . We use the non-dimensional arclength along the ring's undeformed centreline as the Lagrangian coordinate and parametrise the position vector to the ring's undeformed centreline as
We wish to compute the position vector to the deformed ring's centreline.
Here is a sketch of the problem:
Note that we have chosen the Eulerian coordinate axes so that they coincide with the ring's lines of symmetry.
Standard linear stability analysis [see, e.g., G.J. Simitses "An introduction to the elastic stability of structures", Prentice-Hall, (1976)] predicts the ring to become unstable to non-axisymmetric perturbations with an azimuthal wavenumber of at a pressure of
For we therefore expect the ring to deform into a shape similar to the one indicated by the dashed line in the above sketch. As the ring buckles non-axisymmetrically, the material point on the vertical symmetry line (i.e. the material point with Lagrangian coordinate ) moves radially inwards. This makes it an excellent control point for this problem as we can "sweep" through the entire range of the ring's post-buckling deformation by varying its – coordinate from (corresponding to the undeformed, axisymmetric configuration) to (corresponding to a configuration in which the ring is collapsed to the point of opposite wall contact). We note that Flaherty, Keller & Rubinow's analysis [SIAM J. Appl. Math 23, 446–455 (1972)], based on an inextensible beam model, predicts opposite wall contact to occur at an external pressure of
To apply displacement control we change the formulation of the problem from
Determine the position vector to the centreline of the deformed ring, , in terms of the Lagragian coordinate , for a given value of the external pressure .
Determine the position vector to the centreline of the deformed ring, , in terms of the Lagragian coordinate , for an external pressure that is such that the vertical coordinate of the control point is equal to , i.e.
where is given.
The figure below shows computed ring shapes for . They were obtained in two phases: During the first phase of the computation we subjected the ring to a load of the form
where is the outer unit normal on the deformed ring and is a small cosinusoidal perturbation to the external pressure which forces the ring to buckle "in the right direction". The undeformed configuration was used as the initial guess for an initial computation with [or, in the case of displacement control, ]. Subsequently we increased [or decreased ] in small steps, using the previously computed solutions for the displacement field [and ] as initial guesses for the unknowns. This procedure was continued until the ring was collapsed up to the point of opposite wall contact. During the second phase, we set and reversed the procedure to re-trace the deformation to the axisymmetric state.
The figure below illustrates the load/displacement characteristics computed by this procedure. The graph shows the radii of two material points on the ring: The green line shows the radius of the control point; the red line the radius of the material point located at , i.e. the material point located on the ring's second line of symmetry. (Because of the symmetry of the buckling pattern this line may also be interpreted as the load-displacement curve for the control point when the ring buckles in the "opposite" direction). The dash-dotted blue line shows the load/displacement curve when the ring deforms axisymmetrically. In this mode, the radii of the two material points are obviously identical. Finally, the dashed lines shows the load/displacement path when the ring is subjected to a non-zero perturbation pressure, , which deforms it non-axisymmetrically so that even for .
The diagram clearly illustrates the enormous change in stiffness when the ring changes from the axisymmetric to the non-axisymmetric regime. The non-axisymmetric regime emanates from the axisymmetric state (via a supercritical bifurcation at a pressure of , as predicted by the linear stability analysis) and opposite wall contact occurs at , in perfect agreement with Flaherty, Keller & Rubinow's theoretical analysis.
To facilitate the solution of solid mechanics problems with displacement-control techniques,
oomph-lib provides a
DisplacementControlElement, which may be used to add the displacement control constraint to the global system of equations that is solved by Newton's method. Applying displacement control in a solid mechanics problem involves the following straightforward steps:
DisplacementControlElementexpects the load level to be the one-and-only value in the load
Dataobject. We note that computations with a prescribed load are still possible and simply require pinning of the value that represents the load.
Dataobject that stores the load level to the element's external
Datawhose values affect the element's residual vector, but whose values are not determined by the element). External
Datais automatically included in an element's equation numbering procedure. Furthermore, since the elemental Jacobian matrices of
SolidFiniteElementsare generated by finite-differencing, the derivatives of the element's residual vector with respect to the load level are computed automatically. Consequently, the application of displacement control does not require a re-implementation of the solid mechanics elements.
DisplacementControlElementand add it to the
DisplacementControlElementadds the displacement constraint to the global system of equations and thus provides the additional equation required to determine the unknown load level. We note that the
DisplacementControlElementhas two constructors:
Dataobject whose one-and-only value contains the unknown load level. This version of the constructor is appropriate in cases where the load
Datahas already been created elsewhere (as described above).
Dataobject internally and provides access to it via a pointer-based access function.
The driver code discussed below illustrates these steps.
Global_Physical_Variables, used to define the problem parameters and the load function, is very similar to that used in the previous example without displacement control. They main difference is that the adjustable load (the external pressure ) is now defined as a
Data value, rather than a double precision number. This allows it to become an unknown in the problem when displacement control is used.
The main function builds two different versions of the problem, demonstrating the use of displacement control in the cases when the
Data object that contains the adjustable load already exists, or has to be created by the
DisplacementControlElement, respectively. The two versions only differ in the The constructor.
The problem class is very similar to the one used in the previous example without displacement control. Here we store the number of solid mechanics elements in the Problem's private member data since we will add the
DisplacementControlElement to the mesh.
The first half of the constructor is similar to the one used in the previous example without displacement control. We create a
Ellipse with unit half axes) to define the ring's undeformed geometry, and build the 1D Lagrangian mesh. Note that, because of the symmetry of the buckling pattern, we only discretise one quarter of the domain.
Next we apply the symmetry boundary conditions: Zero vertical [horizontal] displacements and infinite [zero] slope at [at ]. (See the previous example for a more detailed discussion of the boundary conditions for beam elements.)
We now distinguish between the cases with and without displacement control:
Dataobject whose one-and-only value stores the adjustable load level (the external pressure), and store the pointer to the newly created
Global_Physical_Variables::Pext_data_pt. Since the value of the external pressure is prescribed (i.e. not an unknown in the problem), we pin the value.
Datadoes not yet exist
Dataobject whose one-and-only value stores the adjustable load. We obtain a pointer to the newly-created
Dataobject from the access function
DisplacementControlElement::displacement_control_load_pt()and store it in
Global_Physical_Variables::Pext_data_ptto make it accessible to the load function
Dataobject that specifies the adjustable load level might already have been created elsewhere in the
Problem. For such cases,
oomph-libprovides an alternative constructor whose argument list includes the pointer to the already-existing load
Dataobject. To demonstrate its use we create a suitable
Dataobject and store it in the
Data(see "Internal" and "external" Data in elements and "global" Data in Problems for a more detailed discussion of this step) and pass it to the constructor:
DisplacementControlElementto the mesh to ensure that it is included in the
Next, we execute the usual loop over the elements to pass the pointers to the problem's non-dimensional parameters the elements. If displacement control is used, we also pass the pointer to the load
Data object to the elements' external
Data to indicate that the element's residual vectors depends on the unknown load. Finally, we set up the equation numbering scheme.
The post-processing function documents the ring shapes and adds the load/displacement characteristics of the two material points on the ring's symmetry lines to the trace file.
We start by opening a trace file to record the load/displacement characteristics, and output the initial configuration.
Next we set up the increment of the control parameter, choosing the displacement or load increments such that the ring's deformation is increased from the axisymmetric initial state to total collapse with opposite wall contact in 11 steps.
Without displacement-control the Newton method can require a large number of iterations, therefore we increment the maximum number of iterations.
We start the parameter study by increasing the ring's compression (either by increasing the external pressure or by reducing the - coordinate of the control point) with to induce buckling in the desired direction.
Then we reset the perturbation pressure to zero and reduce the ring's collapse by decreasing the external pressure (or by increasing the - coordinate of the control point).
In the section The constructor we encountered two different constructors for the
Dataobject that contains the unknown load value is created by the constructor) is most natural for the problem considered here as the load
Datais created specifically for the purpose of allowing displacement control to be applied. It therefore natural to regard the
DisplacementControlElementas being "in charge of" the load
Data, and storing it in the element's "internal Data". (Recall that once
Datais stored in an element's internal
Data, the element becomes responsible for performing tasks such as equation numbering, timestepping, etc. that must be performed exactly once for each
Datavalue in the
Datahas already been created by another element, implying that the other element is already "in charge" of the
Dataobject and performs the equation numbering, timestepping etc. for its values. In that case, the load
Datais regarded as "external Data" for the
In our example code, we simulated the second scenario by creating the load
Data object before calling the second version of the constructor. While this ensures that the load can be regarded as an unknown in the problem, the
Problem remains "unaware" of the additional unknown, as none of the elements in the
Problem's mesh is "in charge" of it. While this is clearly a somewhat artificial scenario,
oomph-lib provides a mechanism for handling such cases: Adding a
Data object to the
Problem's "global Data" by calling
Problem itself "in charge" of this
displ_controlin the main function to
false) and confirm that a non-zero perturbation pressure, is required to induce the ring's non-axisymmetric collapse This shows that the numerical model is not as sensitive to non-axisymmetric perturbations as the theory suggests – roundoff error alone is not sufficient to initiate non-axisymmetric buckling. Use this version of the code to compute the load-displacement characteristics of the ring in its axisymmetric state, i.e. the dash-dotted blue line in the bifurcation diagram shown at the beginning of this document.
Dataand explain why the code fails with a segmentation fault. Use the
Problem::self_test()function to identify the problem. [Hint: What is the default value for a
Dataobject's global equation numbers? What happens if this default is not overwritten? Why is the default assignment not overwritten?]
A pdf version of this document is available.