Stress‐constrained topology optimization of structures subjected to nonproportional loading

This work considers the topology optimization of hyperelastic structures for maximum stiffness (minimum compliance) subject to constraints on their volume and maximum stress. In contrast to almost all previous works, we subject the structures to nonproportional loading, wherein the maximum stress does not necessarily occur at the final load step. As such, the stress is constrained at each load step. The augmented Lagrangian method is used to formulate the optimization problem with its many constraints. In numerical examples, we investigate different load trajectories for the same terminal load and compare the optimized designs and their performances. The results show the importance of considering the entire load trajectory as the load history significantly influences the optimized designs.

the design is prone to stick to the boundary. To overcome this "sticking" effect one may use filter techniques with Robin boundary conditions, see Wallin et al. 8 or "padding" where the design domain is extended, see Clausen. 9 Finally, due to the complexity of enforcing stress constraints, existing TO formulations are, with a few exceptions (e.g., da Silva et al. 10 ), restricted to linear elastic response.
In contrast to the global stiffness or eigenvalue response functions for which a single value quantifies the response function, stress is a field. In our optimization, this means the stress is constrained at every spatial evaluation point in our discretized structure. A common approach to avoid the large number of local stress constraints is to constrain the (maximum) stress and by approximating the maximum stress via the differentiable p-norm, 2 p-mean function, 11 or Kreisselmeier-Steinhauser function. 12 Other methods for accommodating the large number of stress constraints use the augmented Lagrangian (AL) method, see References 10,11,[13][14][15][16][17][18][19][20][21][22][23][24] In addition to these methods, Amir 25 constrains stress by performing an elastoplastic simulation and constraining the effective plastic strain field to be zero which implies that the peak stress does not exceed the yield stress.
Stress-constrained TO studies that consider nonlinear response are rare. Examples include da Silva et al., 10 who design hyperelastic compliant mechanisms and Ivarsson et al., 26 who use an elastoplastic material model and constrain the plastic work. The latter study also investigates load path design dependence, which is used as inspiration for our research. In Reference 26 the primal problem is path dependent, and as a consequence the sensitivity analysis is complex. In this work, the primal problem is quasistatic elastic which simplifies the primal and sensitivity analyzes.
Most stress-constrained TO formulations only consider the stress in the terminal load state. This constraint is only sufficient for structures subjected to proportional loading; however, significant stress levels may emerge at different times if they are subjected to nonproportional loading, which is an important aspect in the design of engineering structures such as engine components, crankshafts, and hinge brackets.
We consider hyperelastic structures subjected to nonproportional loading wherein the entire load history must be considered during the optimization. Recently, Giraldo-Londoño et al. 16 presented an AL approach for the stress-constrained TO of structures subjected to general dynamic loading. Other related works that consider nonproportional loads are fatigue-based TO, compare Zhang et al. 27 and Suresh et al. 28,29 We focus on the nonproportional quasistatic loading to a given terminal load. Notably the peak stress does not necessarily occur in the terminal load state. This is also an important consideration when designing structures that exhibit large deformations, for example, that buckle. Examples of this are seen in the design of compliant mechanisms or bistable elastic structures, compare Wallin et al. 30 and in da Silva et al., 10 which discussed the importance of using nonlinear analysis for structures that endure large displacements. These design problems require a method for constraining the stress in space and in time. In this work we employ the AL method to enforce these many constraints.

KINEMATICS, EQUILIBRIUM, AND CONSTITUTIVE RESPONSE
In our hyperelastic simulations we distinguish between the deformed equilibrium configuration Ω and the undeformed, stress-free configuration Ω o . For quasistatic conditions, the motion of a material particle initially located at position X ∈ Ω o is described by a smooth mapping, that is, the deformation, such that the current position is located at (X) = X + u(X), where u is the displacement field. The deformation gradient, F = ∇ , describes the local deformation, where ∇ is the gradient operator and the Jacobian, J = det(F) > 0, represents the volumetric change of differential volume elements. The local deformation is also described by the right Cauchy-Green deformation tensor C = F T F and the Green-Lagrange strain tensor E = 1 2 (C − 1). The boundary of the body in the reference configuration Ω o consists of the complementary surfaces Ω ou and Ω ot , over which the prescribed displacement u o and the traction t o are applied. Neglecting body forces and inertia, the weak form of the mechanical equilibrium is stated as: find the smooth u which fulfills u = u o on Ω ou and for all smooth u that fulfill u = 0 on Ω ou . In (1), S is the second Piola-Kirchhoff stress which is related to the Cauchy stress via = J −1 FSF T . The virtual Lagrangian strain E in (1), is defined as The hyperelastic material is modelled using an isotropic compressible neo-Hookean constitutive law such that its strain energy density W 0 , is given by where K 0 and G 0 , in the limit of infinitesimal strain, correspond to the bulk and shear moduli. The second Piola-Kirchhoff stress S = 2 W 0 C and the incremental stiffness tensor D = 4 2 W 0 C 2 that follow from the strain energy function W 0 , are provided in the Appendix for completeness.

FINITE ELEMENT DISCRETIZATION
We assume quasistatic loading and apply the load over N ls steps. For each load step, n, the weak form, (1), is discretized using the Galerkin finite element method to obtain the residual equation where a (n) and F (n) ext are the load step n nodal displacement and external force vectors, and r(a (n) ) and F int (a (n) ) the residual and internal force vectors. Dropping the load step superscripts in the remainder of this section for conciseness, the nodal displacement vector is partitioned as , where a p contains the known prescribed displacements on the boundary Ω ou , and a f contains the free displacements, to be determined. The internal force vector is computed by assembling the element Ω e o internal force vectors, as where a e is the element nodal displacement vector, B is the incremental (Lagrange) strain-displacement matrix, compare Crisfield, 31 ⅀ is the finite element assembly operator and N elm is the number of elements in our discretized structure. The internal force vector is similarly partitioned as , where F int,p contains the reaction forces of the prescribed nodes and F int,f contains the internal forces of the free nodes. We similarly partition the external force vector as F ext = [ F ext,p , F ext,f ] and the residual r as r = [ r p , r f ] . To solve the nonlinear residual equation (4), we employ the Newton-Raphson scheme wherein the residual equation r f (a + da) = 0 is linearized about a solution guess a to obtain the unknown displacement increments da f , that solves In accordance with the a = ( a p , a f ) partitioning, the above tangent stiffness matrix K = r a is partitioned as For details of the explicit form of K see for example, Crisfield. 31

MATERIAL INTERPOLATION
In our density-based TO the volume fraction field, is the design field that we optimize. It is discretized to be piecewise uniform over the finite elements, that is, we assign a design variable e , e = 1, … , N elm , for each element in the mesh. We use restriction to formulate a well-posed TO problem. To do so, is filtered to define the filtered design field̃via the partial differential equation (PDE) filter. This filter was originally introduced by Lazarov and Sigmund 7 and later modified by Wallin et al. 8 to mitigate undesired boundary effects. The filtered design field is obtained from the solution to where n o is the outward unit normal vector to the boundary Ω o , Δ is the material Laplacian operator over Ω o and l o and l s are the filter bulk and surface length scales, respectively. Details on how to chose l s and l o are provided in the Numerical examples section. The filter Equations (8) and (9) are discretized using finite elements, see Lazarov and Sigmund 7 and Wallin et al., 8 to obtain the linear system which we solve for the filtered nodal volume fraction vector̃for the given element design variable volume fraction vector .
We employ Heaviside thresholding, compare Guest et al. 32 and Wang et al., 33 to push the filtered volume fraction field̃∈ [0, 1] to the physical volume fraction field ∈ {0, 1} which controls the constitutive behavior. To do this we define = H(̃) where H is the Heaviside step function. But to ensure differentiability for the optimization problem we approximate H with the smoothed Heaviside step function H such that Note that H satisfies H = lim H →∞ H so in fact ∈ [0, 1], that is, regions with intermediate volume fraction where ∈ (0, 1) still exist. During the optimization we apply continuation by increasing H ; we keep = 0.5 constant.
To further reduce the size of intermediate volume fraction regions, we interpolate the material properties via the SIMP scheme, Bendsøe, 34 wherein the strain energy W 0 in Equation (3) is replaced with where Here, p > 1, is the penalization exponent and o is a small residual value which ensures that the mechanical boundary value problem is well-posed. It should be noted that (1) = 1 represents solid material and that in our ersatz material, "void space" is replaced by a compliant material since (0) = o .
As a consequence of the SIMP interpolation scheme, some regions will have a very low stiffness which may cause problems if large deformations occur as highly distorted elements will appear. Wang et al. 33 address this issue by modifying the strain energy whereas Dalklint et al. 35 remove the void elements. In our numerical examples, the deformations are moderate and therefore no special action is required.

OPTIMIZATION FORMULATION
The optimization problem we consider maximizes the stiffness at the terminal load step while constraining the structure's mass and the von Mises equivalent stress of the Cauchy stress (which we will hereinafter call the von Mises stress) throughout the load history, that is, we solve where f is the stiffness objective function, g 0 is the volume and g (n) i are the stress constraints functions defined below.

Objective function
For large deformation, several possibilities to define the stiffness exist. One possibility is the tangent stiffness, see Wallin et al., 36 and another possibility is the end compliance; the latter is used in this work. The end compliance at the terminal load step N ls is unconventionally defined as This is because we use displacement control, that is, the displacement a (N ls ) p is prescribed and therefore we maximize the stiffness by maximizing the reaction force F (N ls ) int,p in accordance with Niu et al. 37 Equivalently we minimize the inverse of (15) 1 We prescribe zero displacement over the surface Ω ou0 and nonzero displacement over disjoint boundary subregions, Ω ou = Ω ou1 ∪ Ω ou2 ∪ … ∪ Ω ouB , and subsequently denote these subregions as Ω oui as p i , i = 0, 1, … , B. To avoid the possibility of vanishing reaction loads over a subregion p i we follow Ivarsson et al. 26 and minimize where the normalization factor C * is chosen so that f = 40, for the initial design. This choice ensures good scaling for the MMA-solver, see Svanberg. 38

Volume constraint
The volume is constrained via where |Ω o | = ∫ Ω o dV is the volume of the design domain in the undeformed configuration and 0 < V * < 1 is the prescribed maximum volume fraction limit.

Stress constraints
To mitigate the singularity problem, we use the relaxed second Piola-Kirchhoff stressS = 2̃( ) W 0 C in the stress evaluation, that is, we use (12) but we replace bỹdefined such that where q < p, compare Bruggi et al. 5 In this way, the stress in the intermediate regions where ∈ (0, 1) is artificially high.
HavingS we compute the relaxed Cauchy stress̃= J −1 FSF T and then the relaxed von Mises stress where [⋅] dev denotes the deviatoric component of the second-order tensor. For later purposes, we also define the relaxed Kirchhoff stress̃= J −1̃a nd̃e ff = J −1̃e ff .
In our stress constrained TO,̃e ff should be less than the threshold stress TH at every point in the structure and for every load step. In our discretized problem, the stress constraint is enforced at each load step n = 1, 2, … , N ls denoted by the superscript, and at each quadrature point in the mesh i = 1, 2, … , N gp denoted by the subscript. The N ls ⋅ N gp stress constraints are formulated as To accommodate the large number of stress constraints, the AL method is employed in this work.

AL method
The AL method replaces the original constrained optimization problem (14) with a sequence of unconstrained optimization subproblems. * The objective function f is replaced by the AL function L, which combines the original objective function with the constraints using quadratic penalty terms.
L( , a (1) , a (2 ) , … , a (N ls ) ; (1 ) , (2) where ⟨⋅⟩ is the Macaulay bracket, vol and the are penalty parameters and vol and (n) i the Lagrange multipliers that enforce the volume and the stress constraints, respectively. A series of optimization subproblems k are solved for the prescribed penalty parameters k Based on da Silva et al. 13 each subproblem k is deemed converged after 50 design iterations, as solving each subproblem strictly up to a prescribed tolerance is computationally inefficient. The penalty parameters are incremented after the completion of each subproblem, k as and the Lagrange multipliers, vol and (n) i are updated after every fifth design iteration † as The detailed algorithm is presented as pseudocode in Algorithm 1.

SENSITIVITY ANALYSIS
The optimization problem (23) is solved using the MMA scheme, Svanberg, 38 which requires gradients of the objective and constraint functions. We solve an unconstrained (with the exception of box constraints) optimization problem (23), so it is not necessary to use an optimization solver that can accommodate constraints. Any first-order optimization algorithm could be used, for example, the steepest descent method with move limits and projection is used by da Silva et al. 13 To evaluate the gradient L we use the adjoint sensitivity scheme rather than the direct method since the number of design variables greatly outnumbers the single Lagrangian function. To derive the sensitivities we define the augmented AL function asL =L( , a (1) , a (2) , … , a (N ls ) ) = L( , a (1) , a (2) , … , a (N ls ) ) + where it is understood that = (̃( )), and a (n) = a (n) ( ). SinceL depends on the displacement in all load steps, n = 1, … , N ls , we need one adjoint vector, (n) f for each load step. DifferentiatingL with respect to the design variables and using the chain rule yields where To annihilate the implicit sensitivity a (n) f , we assign (n) f so that the term in brackets in (29) vanishes, that is, so that where we used the symmetry of K ff . In this way, the sensitivity reduces to The derivativẽis trivially computed from (11) and the derivativẽis annihilated by a second adjoint analysis, compare Lazarov and Sigmund. 7 The sensitivity of the volume constraint is straight forward to compute. It should be noted that although we consider multiple equilibrium states, the adjoint problems (30) are decoupled, that is, we do not have a terminal-value adjoint problem. Thus at each load step n we evaluate a (n) , update L, evaluate (n) and update L .
We now discuss the various derivatives in (31). The external force is zero in our simulations and thus where N e contains the element shape functions used for the interpolation of̃and ′ is the derivative of the SIMP-function. The derivative of the Lagrangian (22) with respect to the displacement vector, that is, the adjoint load in (31) is where The Lagrangian derivative with respect to which appears in the sensitivity (31) is The partial derivative of the compliance f in (17) with respect to a (N ls ) f which appears in (33) is where is the partition of the tangent stiffness matrix corresponding to the free degrees of freedom, and the prescribed degrees of freedom on boundary p i , compare (7). The other compliance partial derivative which appears in (35) is The r T p a p computation is obtained by replacing Λ with a p in (32). The derivatives of the stress functions g (n) i ,̃( eff are also required in (33) and (35). Using the chain rule on (20) results iñ where we use C (n) a (n) f = 2B (n) , and define the stress like quantity where 3̃e ff and C = J −2∕3 C is the isochoric Green-Lagrange tensor. The derivation of S is provided in Wallin et al. 30 A straight forward derivation reveals the partial derivativẽ

NUMERICAL EXAMPLES
The presented formulation is used to design three maximally stiff structures subject to volume and maximum von Mises stress constraints. We restrict the examples to two-dimensional problems as the aim of this work is to investigate the optimization algorithm and the impact of nonproportional loading. Extending the framework to three-dimensional comes with an additional computational cost but is straightforward. All problems are solved for a maximum volume fraction V * = 0.35, initiated with a uniform design field = 0.35 and solved via the MMA with standard parameter settings by Svanberg, 39 except for an explicit control of the moving asymptotes. We restrict the maximum distance between an asymptote and the design variable to 0.1, as discussed in Verbart et al. 40 and for the last 100 iterations we further restrict the move limit, that is, the maximum change for a design variable to be 0.05. The material penalty parameters are held constant during the optimization at p = 3 and q = 0.5. A continuation scheme is applied to the Heaviside threshold parameter H such that we initiate H to 1 and when a "clear structure" has emerged, that is, after 250 iterations, it is increased by H ← 1.06 ⋅ H , every 10th design update until the final value H = 12 is reached. The initial values of the penalty parameters are set to 0 = 5000 N ls N gp and 0 vol = 10 where N ls N gp is the number of stress constraints. The Lagrange multipliers are initiated to zero. The penalty parameters are increased by continuation according to (24) until the maximum values, max = 1 and max vol = 300, are reached. The optimization terminates 100 iterations after the final value of H is reached since no significant design changes are seen thereafter.
The structures are discretized using four-node bilinear quadrilateral plane strain elements with the standard four-point Gauss-Legendre quadrature scheme. The bulk filter length scale is l o = r∕(2∕ where l e is the element side length. Based on the work of Wallin et al., 8 we equate l s = 0 over the boundary regions p 0 , p 1 , … , p B and boundary locations within a distance 2r from these. We equate l s = l o over the remainder of the boundary. Furthermore, due to the applied displacement boundary conditions, we increase the effective stress threshold over regions extending a distance r from p 0 , p 1 , … , p B , to be 50% larger than for the rest of the structures.
The isotropic material is modeled via the Young's modulus E = 100 GPa and Poisson's ratio = 0.3 which provides the bulk and shear moduli;

L-bracket
The first problem we consider is the bench mark L-bracket example, compare Figure 1. The structure with L = 35 m is subjected to proportional loading by prescribing the displacement,ū y = 5 m, over the upper portion, p 1 , of the right surface, and zero displacement over the top face p 0 . We limit the effective von Mises Cauchy stress to TH = 2.4 GPa. The design domain is discretized with 19,600 quadratic elements with a side length of l e = 0.2 m. As seen in Figure 2, the optimized Figure 2A design avoids the stress concentration at the reentrant corner. The peak stress occurs over the p 1 loading surface which we attribute to the boundary condition. For this reason we are more interested in the peak stress eff = 2.41 GPa in the vicinity of the reentrant corner and the fact that much of the structure is fully stressed, see the box in Figure 2B. A similar design was conceived in Ivarsson et al. 26 by constraining plastic work.
The optimization convergence history for the L-bracket is shown in Figure 3. The effect of the continuation strategies of the Lagrange multipliers and penalty parameters are clearly seen in Figure 3A. One can also notice the initiation of the Heaviside threshold parameter continuation at design iteration 250 in Figure 3C as the peak stress starts to oscillate. This is because the peak stress is located on the smeared boundary region, which is highly effected by the Heaviside projection.

Double L-bracket
In the next example we consider the double L-bracket, shown in Figure 4 with the side length L = 150 m. We prescribe the displacements, u y1 and u y2 , on the center of the left p 1 and right p 2 boundaries and clamp the top surface p 0 . Two load patterns are investigated to illustrate the effect of nonproportional loading. The loads are applied in eight steps as shown in Figure 5. For the proportional load pattern Figure 5A, we simultaneously increase the displacement on both sides until the final load is applied. In the second, nonproportional load pattern Figure 5B, we first apply the displacement over the right boundary until it reaches 75% of the terminal load, and then apply the displacement on the left side. Both loadings have the same terminal state, that is, u y1 = u y2 = 17.6 m =ū. The effective von Mises stress limit is TH = 2.3 GPa. The design domain is discretized with 19,800 quadratic elements with a side length of l e = 1 m. Designs optimized for the two load paths are illustrated in Figure 6. The most pronounced difference appears in the center region. For the proportional load design we obtain a sharp V-shaped region whereas the nonproportional design contains a smooth U-shaped region to avoid high shear stresses. These designs are also similar to those produced by Ivarsson et al. 26 The von Mises stress distributions at the terminal load step, shown in Figure 7 indicate that stress concentrations are avoided. The peak stress of the proportional and nonproportional load designs violate the stress constraint by 0.8%, and 0.5% and their end compliances are 46.63 and 45.78 GNm. The design optimized for the proportional loading is approximately 2% stiffer under the final load step than if optimized for the nonproportional loading even though the terminal loads are the same and we only consider the end compliance in the optimization. This difference is due to the stress constraints. Figure 8 shows the peak stress evolution for the nonproportionally loaded design after each load step corresponding to design iterations 10 (A), 200 (B), and 780 (C). We see that the peak stress for the last five load steps approaches the threshold stress in the converged design. The stiffness of the nonproportionally loaded design presumably suffers in order to satisfy these additional stress constraints. Next, the nonproportional load is applied to the proportionally loaded design, and it is observed that the maximum stress level occurs in the fourth load step. Figure 9 compares the von Mises stress distribution after the fourth load step for the two designs. As expected we obtain significant stress concentrations for the proportionally loaded design. In fact the threshold stress is exceeded for a majority of the load steps, compare Figure 10. The boxes in Figures 7 and 9 show the stress distributions in both designs at load steps 4 and 8. Note that the location of the maximum stress for the proportionally loaded design subjected to the nonproportional loading is at the V-corner.

Cross beam
In the final example we analyze a cross structure with a side length L = 150 m, compare Figure 11. We limit the effective von Mises stress to TH = 4.4 GPa.

F I G U R E 11
Geometry and boundary conditions for the cross beam.

F I G U R E 12
Optimized designs, for the load patterns #1, #2, and #3.
We prescribe displacement on the left, p 1 , right, p 2 , and bottom, p 3 , middle boundaries as shown in Figure 11 and clamp the top p 0 surface. To demonstrate the load path influence we vary the boundary conditions in three ways, compare Figure 12. For the #1 proportional case we apply the loads on all boundaries simultaneously, that is, u y1 , u y2 , and u x3 are increased at the same rate. In the #2 load path we first simultaneously apply the loads on the bottom, p 3 , and the left, p 1 , side, and then on the right, p 2 , side. In the #3 load path we first displace the bottom, p 3 , boundary, then the right, p 2 , side and finally the left, p 1 , side. The fully loaded state corresponds to prescribed displacements u y1 = u y2 = u x3 = 14.4 m =ū. The design domain is discretized with 25,200 quadratic elements with a side length of l e = 1 m.
The optimized designs for the three load patterns are shown in Figures 12 and 13. The most pronounced differences between the designs is the small shift of the connection to the top boundary and the left side of the structure. It can also be seen that the Figure 13C design is slightly thicker in the bottom region. The von Mises stress distributions for the optimized designs are presented in Figure 14 and the von Mises peak stress and end compliance are presented in Table 1 for the final load step.
As seen in Table 1, the peak stress of the three designs at the terminal load step is similar. The end compliance, that is, the objective function varies considerably more and it is seen that the stiffest design is obtained for the proportional load pattern #1. This is not surprising based on the results of the previous double L-bracket example.
To further demonstrate the possibly catastrophic consequences of ignoring the load trajectory in the optimization, load path #3 is applied to all designs, compare Figure 15 and Table 2. As expected, the maximum stress for the Figure 13C design optimized for #3 load path, barely exceeds the stress limit TH whereas the maximum stress for the Figure 13A,B designs significantly exceeds TH at the several load steps.  Finally, in the boxes in Figure 16 we present the location of the peak stress in the center regions. We note that the maximum stress actually occurs at the loaded boundaries but these stress concentrations are the consequence of the boundary conditions which are unavoidable.

CONCLUSIONS
A TO framework based on finite strain hyperelasticity for maximizing the stiffness while constraining the mass and the maximum von Mises effective stress for the entire load trajectory is presented. The optimization problem is solved using the AL method which has proven to be effective for accommodating a large number of stress constraints and it is regularized by imposing a minimum length scale via the augmented version of the PDE filter to avoid artificial boundary effects. Crisp, that is, black and white designs are promoted by the use of Heaviside thresholding. The adjoint sensitivity scheme is used to derive the sensitivities, and the design updates are obtained from the method of moving asymptotes.
In contrast to TO problems that consider plasticity or dynamics, the state problem is decoupled due to our quasistatic loading assumption. From a computational point of view, this assumption simplifies the sensitivity analysis as we do not need to save or recompute the stiffness matrices. The computational cost is, however, similar as the number of linear equations solved is approximately the same. However, if one were to use a linear elastic model, all load scenarios could be solved effectively in parallel using a direct solver, which would result in a substantial speed up as discussed in Koppen et al., 41 compared to using our hyperelastic model. In this work, several equilibrium states has to be solved separately using Newton's method, which contributes to the computational cost.
To validate the TO designs, da Silva et al. 15 used a body-fitted mesh without void material to compute the response. They concluded that the peak stress levels can be underestimated by up to approximately 12% in density-based TO compared to a body-fitted mesh. This important matter requires further research.
The numerical examples show that the load path significantly influences the optimized design and that to prevent failure, the entire load history should be considered as it is not always known where, or when, the peak stress will appear.
The optimization problem is solved using the in-house Fortran 90 implementation summarized in Algorithm 1 together with the MMA code from Prof. Krister Svanberg. Linear equation systems are solved using Intel MKL Pardiso 2020.

ACKNOWLEDGMENTS
This work was performed under the auspices of the US Department of Energy by Lawrence Livermore Laboratory under Contract DE-AC52-07NA27344. The Swedish energy agency (grant number 48344-1) and eSSENCE: The e-Science Collaboration (grant number 20206:1) are also gratefully acknowledged as well as the Swedish Research Council (grant number 2021-03851). Last but not least, the authors would like to thank Dr. Niklas Ivarsson for helpful discussions.

CONFLICTS OF INTEREST STATEMENT
The authors declare that they have no conflict of interest.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. ENDNOTES * The PDE governing equations are constraint equations. However, we do not incorporate them as constraints in (22)