Stress resultant plasticity for plate bending in the context of roll forming of sheet metal

Modeling of roll forming process of sheet metal requires efficient treatment of plastic deformations of thin shells and plates. We suggest a new general approach toward constructing the governing equations of bending of an elastic‐plastic Kirchhoff plate on the structural mechanics level. The generic function of the isotropic hardening law is formulated in terms of variables, which are meaningful for a material surface, namely resultant bending moments and dissipative work. This function is then identified by comparison of the exact solutions for a uni‐axial or isotropic bending experiment within the structural model and with the continuum model of a through‐the‐thickness element. We validate the approach against both the fully three‐dimensional simulations as well as the results of the traditional elastic‐plastic plate analysis, the latter one featuring the continuum laws of plasticity treated in the integration points along the thickness direction. We considered benchmark problems from the literature as well as a prototype for the roll forming problem. Besides good agreement regarding the shape of plastic zones and force‐deflection characteristics, experiments also demonstrate higher computational efficiency of the new stress resultant model.


INTRODUCTION
In the present paper, we introduce a novel approach to modeling plasticity in thin plate-and shell-like structures at bending, which is computationally more efficient than traditional ones. The ultimate aim of our simulations is set on the roll forming of sheet metal, where the proposed model of thin elastic-plastic plates shall play a fundamental role.
Roll forming belongs to the class of cold forming manufacturing processes. The basic idea is that an initially flat sheet metal strip is transformed to the desired product with a different profile (shape of the cross section), see References 1,2. The strip is continuously fed through the consecutively arranged forming passes. Each forming pass consists of a set of contoured rolls, which incrementally contribute to the finally obtained deformation, thus producing the desired shape mainly by bending. Therefore, roll forming is essentially different to rolling: no significant thickness deformations occur, the membrane strains remain small, and the final product shape is acquired by irreversible plastic bending deformations.
Roll forming is a highly economical manufacturing process for the mass production of high-quality profiles. Their complexity can vary from simple profiles (e.g., U-profiles) to highly sophisticated ones. While being very important for practical applications, predicting the final shape for a given assembly and design of the forming passes is currently an extremely time-consuming computational task, because traditionally continuum models with Lagrangian kinematic formulation of the finite element scheme are put into practice. Because of the axial motion of the deformable structure through a series of contact zones, it is computationally more efficient to use the mixed Eulerian-Lagrangian kinematic description with a nonmaterial finite element formulation as discussed in . Previous attempts to apply this approach for roll forming featured solid continuum elements 2 and were not very successful from the point of view of numerical efficiency. Faster simulations are expected with structural mechanics models of plates and shells, which require a consistent and efficient treatment of elastic-plastic behavior at bending dominated deformations.
The proposed stress-resultant model of plasticity allows to reduce the computational complexity of this nonlinear problem to a great extent.
In the present paper we consider static deformation of plates with the focus on the plasticity formulation. A widely utilized approach to treat plasticity in thin bodies avoiding the full three-dimensional analysis are structural models of beams, plates, and shells, in which the continuum mechanics laws of elastoplasticity are treated in chosen points of a through-the-thickness element. The kinematic assumptions of the structural model relate the local strains with the overall strain measures of the structure (e.g., curvature of a beam or of a plate), and the stress resultants (forces and moments) follow by integration. A good correspondence to full three-dimensional solutions is usually observed for thin structures with sufficient number of integration points, see discussion in Reference 6. A further step toward computational efficiency can be taken, if the constitutive relations directly connect the stress resultants to the local strain measures of the structural theory (e.g., a relation between the bending moments and curvature for a beam or for a plate) and thus the through-the-thickness integration becomes obsolete.
In the literature, one finds several attempts to describe the elastic-plastic behavior of plates and shells completely on the level of structural mechanics, mainly in the view of computational aspects. Practically all formulations make use of the associated flow rule and an extended form of the Ilyushin yield criterion. 7 The differences lie in the underlying structural plate model (with or without shear and membrane deformations), in the details of the finite element kinematic approximation and, what is most important, in the formulation of the isotropic and/or kinematic hardening law. 8 Moreover, the choice of the internal hardening variables, which enter the evolution law, is also not unique. Thus, Simo and Kennedy 9 introduced a hardening potential that is quadratic in terms of the hardening variables, which implies linear hardening behavior. Two hardening laws are obtained by using that potential. One for isotropic hardening, where the internal hardening variable takes the form of a dissipative function, which is defined as a generalization of the equivalent plastic work. The other introduced law considers kinematic hardening and is an extension to the Prager-Ziegler kinematic hardening law of the classical J 2 -flow theory. The practical application of these laws requires, however, material parameters such as isotropic and kinematic hardening moduli of the structural mechanics theory to be known. Ibrahimbegović and Frey 10 considered the Ilyushin yield criterion adjusted by an additive hardening term in the form of a product of a constant isotropic hardening modulus and the effective plastic strain as an internal variable within a kinematically simple model of a Reissner-Mindlin plate.
While this yield criterion is simple and perfectly applicable, it does not account for the spreading of the plasticity through the thickness of the plate and is useful only for the analysis of limit loads. Crisfield 11 suggested an evolution law for an additional state variable, which reflects the growth of the plastic zone in a through-the-thickness element. Dujc and Brank implemented this strategy in Reference 12 and demonstrated satisfactory results for several benchmark problems, one of which we also consider in the present work. The same authors further extend the approach toward more complicated shell kinematics in Reference 13 making use of a yield criterion function that corresponds to the Ilyushin-Shapiro multisurface flow criterion in three dimensions. Both isotropic and kinematic hardening are accounted for in the proposed yield criterion. The internal hardening variables of their model feature the plastic strain, another scalar parameter which controls the isotropic hardening mechanism and strain like parameters which control the kinematic hardening mechanism.
A different approach of considering hardening in terms of hardening moduli (which are often constants) is pursued by Chou et al. 14 The Ilyushin yield criterion is again applied; however, upon initiation of plastic flow the yield criterion function takes the value of a normalized equivalent stress resultant Σ, which is variable because isotropic hardening is considered. Σ is handled as a function of the normalized equivalent generalized plastic strain and is formulated in terms of a power law, where the exponent n is called hardening exponent, which is viewed as a material parameter and must be known a priori in order to apply the structural model.
In the present paper we focus on the Kirchhoff plate model with the associated flow rule and the Ilyushin yield criterion, which is adjusted in order to account for isotropic hardening; see Reference 15 for a discussion of more general plasticity models in application to sheet metal forming. Similar to Reference 14 we propose, that a through-the-thickness element yields, when the yield criterion function equals a limiting value k 2 M , with k M being the maximal equivalent elastic moment. In contrast to many other approaches, we consider the limiting value k M as a function of the dissipative work A p in the given point of the plate, which fully determines the isotropic hardening of the structural model. Furthermore we do not assume any particular form (e.g., power law) for k M (A p ), but rather identify the analytic function by solving simple problems like uni-axial bending. The identification relies upon the continuum model of a through-the-thickness element: equations of the three-dimensional elastoplasticity are solved for the points on a line in the thickness direction under the kinematic assumptions of the Kirchhoff plate theory. Integrating the three-dimensional stresses and dissipative work over the thickness, we finally obtain the sought for hardening function of the structural model. We do not consider kinematic hardening and thus are not able to account for the Bauschinger effect. 8 Because reverse plasticity is usually avoided in the practice of roll forming, the present model shall be sufficient for this field of application. In numerical experiments, we observe good agreement regarding the shape of plastic zones and force-deflection characteristics. Moreover, the computational advantage compared to conventional plasticity models is pronounced.

GENERAL MODEL OF AN ELASTIC-PLASTIC KIRCHHOFF PLATE
A plate is considered as a material surface in the framework of the so-called direct approach. 16,17 The absence of the thickness coordinate makes it very easy to arrive at the complete theory. Formulating the constitutive conditions requires, however, identification by comparison to three-dimensional solutions, which we will accomplish in the subsequent sections. The single strain measure for a linear Kirchhoff plate is the tensor of curvatures with ∇ being an in-plane differential operator and w being the transverse deflection. With Cartesian coordinates x, y in the plane of the plate and corresponding basis vectors e x , e y the invariant differential operator takes the form such that the symmetric in-plane curvature tensor x y (e x e y + e y e x ) + 2 w y 2 e y e y = x e x e x + xy (e x e y + e y e x ) + y e x e x , is defined by three components. In the following we use the direct tensor notation with a dyadic product implied when two vectors are written next to each other, thus omitting the symbol ⊗, which is sometimes used in the literature. This utilized variant of the mathematical notation is common in the theory of elasticity and structural mechanics, see for example, Reference 18. The elastic law for the isotropic plate is given by in which M is the plane tensor of moments acting in the plate, I is the plane identity tensor, is the Poisson ratio and D is the flexural stiffness. Note that M shall be considered as one of the stress resultants in the three-dimensional body of a plate. The tensor of elastic curvatures e is identical to the total curvatures as long as the deformation is purely elastic. Using the principle of virtual work, one obtains the equilibrium conditions in the form Here Q is the vector of transverse forces and q is the external load per unit surface of the plate, which acts in the transverse direction. The equilibrium equations and the corresponding boundary conditions for M and Q are equivalent to the variational principle of minimum of total potential energy of the plate, which we use below for the finite element analysis.
We account for inelastic deformations by means of the additive decomposition of the total geometric curvature into the elastic and the plastic parts, An evolution law for p is required to close the system of equations. Additional relations will need to be introduced for this sake, namely: yield criterion, hardening law, and flow rule. Each of them is formulated in stress resultants which leads to our general model of an elastic-plastic Kirchhoff plate. Furthermore we focus on a plate model with isotropic hardening. Not taking the phenomenon of kinematic hardening into consideration, we will not be able to account for the Bauschinger effect, 8 which is important in case of reverse plasticity. Nevertheless, we decide that the consideration of isotropic hardening is sufficient for our purposes in the field of roll forming, as typically in roll forming, the sheet metal is bend progressively under several, subsequently acting, shaped roll pairs and thus no plastic reverse bending is expected.
We start by introducing the yield criterion function f (M) and its derivatives. Motivated by the consideration, that the plane stress tensor in the three-dimensional body of the plate is proportional to M and making use of the von Mises yield condition, we formulate the yield criterion function as an invariant quadratic form the double dot multiplication "⋅⋅" has a meaning of a full contraction of two tensors, see which corresponds to the definition given by Ilyushin. 7 The above expression differs from the conventional yield criterion function in the three-dimensional plasticity at plane stress 9,10,19 only by substituting the moment components instead of the components of the stress tensor. We note that (7) In each point, the plate deforms elastically as long as f does not exceed certain threshold: f < k 2 M . This boundary is determined by the maximal equivalent moment k M , at which plastic flow is initiated such thaṫp ≠ 0. While the ideal plasticity model with constant yield surface k M = const is perfectly valid in the three-dimensional case, it is inappropriate for beams and plates because the actual plastic zones in a through-the-thickness element will grow as bending progresses, as discussed in the subsequent sections. We describe this phenomenon by means of the isotropic hardening model: the yield surface expands as plastic flow takes place and thus k M is a function of an internal hardening variable. The latter we choose as the dissipative work A p : k M = k M (A p ). The time rate of the dissipative work done in a given point of the plate equals the dissipative powerȦ Determining the specific form of the function k M (A p ) will be the main goal of the subsequent section. The evolution law for the plastic curvature p , which is needed to complete the general formulation of the elastic-plastic model, is determined by the postulate of maximum plastic dissipation 8 and is formulated as an associated flow rule in the forṁp such that the flow direction is normal to the yield surface. The nonnegative time derivative of (which is sometimes called the plastic consistency parameter) follows from Equations (7) to (11). First we reformulate (10) by using (11) and (7): During the active plastic flow process we remain on the yield surface, and the following relations are valid: Substituting (7) and (12) into (13) we explicitly finḋ: The expression is valid only if it is positive, otherwise the active plastic flow ends immediately and the plate keeps deforming elastically in this point with the reached plastic curvature p .
Assuming the time history of the moment tensor M(t) in a given point of the plate to be given, we can now determine the kinematic response of the plate by taking the following steps: 1. We find the plastic work done by dissipative forces A p (t) by solving the first-order differential equation, which results after substitutinġfrom (14) into (12). In general, one needs to take the unloading condition into account:̇= 0 as soon as the right-hand side of (14) becomes negative. Further plastic flow is initiated after an elastic phase as soon as the first equality in (13) is again fulfilled and the yield surface is reached. 2. Now thaṫ(t) is known, we integrate the flow rule (11) over time.
In a practical computation, neither M(t) nor the kinematics of deformation are known in advance. Moreover, the problem becomes nonlocal as the moment distribution over the plate at static equilibrium depends both on the external loads as well as on the distribution of plastic curvatures. The time integration scheme applicable in the framework of a finite element simulation will be discussed below in Sect. 5.

IDENTIFICATION OF THE ISOTROPIC HARDENING FUNCTION BY COMPARISON TO THREE-DIMENSIONAL SOLUTIONS
We complete the theory and identify the isotropic hardening function k M (A p ) by performing simple thought experiments for a continuum model of a through-the-thickness element of a plate with thickness h, which features integrated three-dimensional equations of elastoplasticity. At doing so, we tacitly assume that integral quantities, such as curvatures, bending moments and dissipative work are equal in both, the structural and the continuum models.
We begin with the case of the uni-axial bending with the moment Now the plate behaves similar to a beam at pure bending and responds with the curvature = x e x e x + y e y e y ; the negative curvature in the transverse direction y appears already during the elastic phase because of the Poisson effect.
Our aim is to obtain the moment-curvature characteristic with the structural model and using the continuum solution to identify k M (A p ) by comparison. We assume the axial bending moment M(t) to be monotonously increasing in time, such that no reverse bending occurs, and begin with the continuum solution for the three-dimensional body of the plate under the assumptions of Kirchhoff's theory regarding the kinematics of deformation and orders of magnitude of various stress components, see Reference 20 for the asymptotic justification. While the simple theory of elastic-plastic bending of a beam can be found in the literature (see e.g. Öchsner 21 ), here we briefly recall the key points, needed for the subsequent analysis.
The conventional Kirchhoff plate kinematics result into a linear distribution of the axial strain component x over the thickness: here z is the coordinate in the thickness direction with the origin in the middle of the plate. This linearity is the asymptotic consequence of the three-dimensional equations of compatibility of strains 20 and thus holds independently from the material behavior. In the elastic range this also means a linear distribution of the stress tensor with E being Young's modulus. All components of the stress tensor but the axial one vanish because of the plane stress assumption of the plate theory and because the plate is free to bend in the transverse direction. Simple integration relates the bending moment to the stress distribution: Computing the integral in the elastic range, we find the elastic relation between the moment and the curvature, which also provides us with the flexural stiffness coefficient in (4): Plastic flow is initiated at the upper and lower surfaces as soon as the stress reaches the yield stress k by the absolute value, from which we find the maximal elastic curvature and the maximal elastic moment to be * = Now we let the bending moment grow beyond the elastic limit, M(t) > M e * , when two plastic zones −h∕2 < z < −z p and z p < z < h∕2 begin to develop and the elastic core −z p < z < z p starts shrinking, as the boundary between them progresses from h∕2 at = * toward 0 at → ∞. The expression for z p follows from the conditions, that the elastic law (18) holds in the elastic core and that | (z p )| = k, as plasticity is just being initiated in this fiber. The distribution of (z) in the plastic zones follows from the corresponding material characteristic. The analytic derivation below makes use of the assumption of the ideal elastic-plastic continuum material behavior, we discuss the implications of more realistic stress-strain curves in the end of this section. Under the stated assumption, the stress equals the yield stress in the plastic zones, = ±k, and the total bending moment in the plastic regime is computed according to (19) as we used the symmetry of (z) and substituted z p from (23). The moment-curvature relationship is thus defined piecewise, and it has a horizontal asymptote While the traditional consideration of the elastic-plastic beam bending usually stops at this point, we keep on computing and find an inverse relation to (25) for the total curvature as a function of the bending moment: The total curvature = e + p is a sum of the elastic part e = 12M∕(Eh 3 ) and of the plastic curvature p , see (6) and (4). Further analysis requires p being expressed as a function of the total curvature: Vanishing of p in the elastic range is obvious, and the expression for ≥ * follows after subtracting e from both sides of (27) and further substituting M = M p ( ) from (24). Now we proceed to compute the dissipative work per unit area of the plate A p done in the plastic zones. The dissipative power of the internal forces is the product of the stress tensor with the time rate of the plastic part of the strains. This power equals the time rate of the dissipative work A p 3 per unit volume: The computation is simple in the present case of ideal elastoplasticity, as the stress remains constant. For certainty we consider the upper plastic zone z p < z < h∕2 with compression, = −k. At known stress and total strain it is easy to find the plastic part using the elastic relation: and with (17) the total dissipative work in a through-the-thickness element follows to The expression makes only sense in the elastic-plastic region ≥ * . The function A p ( ) vanishes together with its first derivative at = * and grows linearly at → ∞. The inverse relation (A p ) is then found to be Now we switch out attention to the structural plate model from Section 2 and evaluate it for the problem of uni-axial bending at hand. Using the same Cartesian coordinate system as before, we write the elasticity law (4) for the bending moment (15), invert it and find the total curvature in the axial direction to the other components are of little interest. As already stated above, the yield criterion function (8) takes now the simple form: The governing differential Equations (11), (12), and (14) result intȯ The component p xy vanishes. We furthermore eliminatėp y and solve foṙp x andȦ p , thus obtaining expressions, which correspond to our expectations in the simple case of beam bending: The denominator in both time derivatives is the same and contains the isotropic hardening function, which we seek to identify. The expression in the denominator can be written more conveniently by introducing a new function which allows us to reformulate (36):̇p DividingȦ p bẏp x , we eliminate the explicit time derivatives and arrive at the differential relation, which will hold in the structural model regardless of the choice of the functions k M and : the relation makes sense only in the elastic-plastic range. Now the question arises whether this relation holds for the continuum model of a through-the-thickness element. Computing the derivatives of A p and p with respect to according to (31) and (28), we find which differs from the expression (24) for the bending moment in the plastic range M p ( ). Although this observation means that the exact correspondence of both models cannot be reached, we nevertheless continue the process of identification, such that the most important characteristic of the continuum model, namely the moment-curvature relation, is described by the structural model as accurately as possible. To this end we express the unknown function from the first equality in (38): We evaluate the right-hand side in the elastic-plastic regime of the continuum model to according to (24), (28) and substitute = (A p ) from (32), which finally results into the sought for expression Integrating (43) with respect to A p and taking into account the initial yield condition k 2 M (0) = (M e * ) 2 at the beginning of plastic deformation with A p = 0, we obtain the yield criterion function to be The obtained functions and k M determine the evolution of the yield surface and thus the growth of the values of the yield criterion function f , see (8), (13), this corresponds to the isotropic hardening effect. As discussed earlier, the continuum and the structural mechanics models cannot be identical, and the accuracy of the approximation depends on the question, how well the differential equation for A p is fulfilled. Before we answer this question in the subsequent section with the help of simple numerical simulations, we shortly discuss the second option to identify the hardening behavior based on the thought experiment of isotropic plate bending.
We follow the same methodology as above, but the tensor of bending moments and the curvature tensor take now the form the same variables obtain a new meaning now. The plane parts of the stress and the strain tensors in the continuum model are isotropic, In the elastic range holds The yield condition at general plane stress is known to be 7 which in our case x = y = , xy = 0 reduces to the same relation as above: The maximal elastic curvature and the height of the boundary between the elastic core and the plastic zone, after plasticity is initiated, follow to * = The total bending moment again follows as (19). Knowing the stress distribution in the elastic core (47) and in the plastic zones (49), we compute the integral and find The maximal elastic moment M e * and the plastic one M p * are the same as above in the uni-axial case. The plastic curvature as a function of the total one p ( ) follows by inverting (51) and subtracting the elastic To compute the dissipative work in a through-the-thickness element A p , we again begin with the dissipative work per unit volume A p 3 . According to the first equality in (29), it takes now the form because both the stress tensor and the tensor of plastic strains have two identical components in the plane of the plate, and = −k as we again focus on the upper plastic zone with uniform compression. The plastic strain follows from the elastic relation (47), in which we now replace the total strain by its elastic part: Acting further as when deriving (31) and (32), we first find A p ( ) by integration over the upper plastic zone z p ≤ z ≤ h∕2, and then invert the relation to obtain the curvature as a function of the dissipative work (A p ).
Proceeding, we treat the structural plate model in a full analogy to the above derivation of (33)-(36). The yield criterion function f retains its form as in the uni-axial case, the differential relations for the plastic curvature and the dissipative work look noẇp Again, either the first relation or the second one may be made compatible with the continuum model by the appropriate choice of (A p ). We choose the differential relation foṙp and express The analytic integration with respect to A p is again possible. Accounting for the initial yield condition, we find a lengthy expression for the square of the hardening function k 2 M (A p ). It is natural to expect, that the accuracy of solutions of the structural plate model with each of the two hardening functions, identified, respectively, for the cases of uni-axial and isotropic plate bending, shall depend on whether the actual bending state is closer to the uni-axial or to the isotropic one.
Finally, we emphasize the fact, that the previously followed procedure to identify k M (A p ) is not restricted to ideal elastic-plastic continuum. Generally, empirically estimated isotropic hardening laws (see Reference 22) find use in roll forming simulations (thus, the isotropic hardening law of Swift is applied in Reference 23). Hardening behavior of the stress-strain characteristic of the material of the plate can also be incorporated in our model by virtue of certain modifications of the above procedure, which we now consider for the case of uni-axial bending. The stress tensor again takes the simple form of (18). Subsequently we introduce an isotropic hardening law, for example a power law according to Swift: where n is a material parameter named strain hardening exponent, which is physically limited to only assume values in the range of 0 ≤ n ≤ 1. The limiting case n = 0 corresponds to the ideal elastic-plastic material behavior, whereas n = 1 corresponds to purely elastic behavior. The moment-curvature relation of the continuum model follows with the same procedure as before, simply by integration of the stress distribution along the thickness. In the general relation (25) we thus find a new expression for the moment in the plastic range: The expressions for the maximal elastic curvature * and for the boundary between elastic and plastic region z p are identical to the uni-axial case without hardening. Furthermore it is easy to verify that for the case of n = 0 all other expressions will also be identical to the expressions above. At n > 0, however, closed form expressions for the sought for functions k M (A p ) and (A p ) are no longer available, because the moment-curvature relation (57) is not solvable analytically. Nevertheless it is still possible to find the plastic curvature and the dissipative work as functions of the total curvature: Equations of the general structural plate model remain unchanged. We again identify the function according to (41), which is available in dependence on the total curvature: ) .
As our implementation rests upon the dissipative work A p being the internal hardening variable, we deal with this problem by first numerically solving the Equation (59) for and then substituting this numerically obtained characteristic in (60), such that (A p ) can be evaluated. The necessary isotropic hardening function k M (A p ) follows by integration, which is now only numerically possible. However, as (A p ) is singular in the vicinity of the initial state A p = 0, numerical integration proves to be challenging. We approach this issue by introducing the function see (37). Subsequently, by the use of (59) and (60) it is possible to analytically derive the lengthy expression for the derivative (61). The thus gained function is valid for ≥ * and is no longer singular. Therefore, numerical integration can be carried out very easily and provides us with k M ( ). The sought for dependence k M (A p ) follows with the use of the numerical characteristic (A p ).

ELASTIC-PLASTIC RESPONSE OF A THROUGH-THE-THICKNESS ELEMENT
Prior to proceeding to finite-element simulations for a whole plate, we validate the derived local constitutive relations by comparing curvature-moment characteristics of a through-the-thickness element, computed for a particular parameter set using both the structural mechanics equations as well as continuum mechanics solutions. The latter continuum mechanics solutions are available analytically for simple cases of uni-axial and isotropic bending, but they result from numerical integration of systems of differential equations when the direction of the bending moment is changing in time. The used set of parameters for all subsequent comparisons is provided in Table 1; SI system of units will be used throughout the paper.
We begin with a uni-axial bending experiment and use the previously derived continuum solution of (25). A solution for the structural plate model follows after the integration of the system of ordinary differential Equations (36).
A moment-driven simulation is considered, such that the bending moment M(t) is known and monotonously increasing in time. Remaining with the assumption of ideal elastic-plastic continuum, we test both previously identified hardening functions, derived respectively for uni-axial and isotropic bending. Subsequently we solve the system of differential equations numerically and obtain p x (t) and A p (t). Thus, by inverting the relation, and by considering the elastic range by means of the elastic law, given in (4), we acquire the complete moment-curvature characteristic. Finally, all three moment-curvature characteristics are depicted in Figure 1. Besides the general good correspondence of both structural models to the continuum one, we also note that the solution with the hardening function, identified for the uni-axial bending problem, expectedly shows slightly better correspondence, in particular in the beginning of the plastic stage.
We also perform a similar experiment for the case of isotropic bending. The continuum solution for the moment-curvature characteristic is given by (51). And for the structural plate solution we need to solve the previously introduced system of ordinary differential equations (54) under consideration of (37). The three moment-curvature characteristics are plotted in Figure 2.
We again observe, that the solution with the corresponding hardening function, obtained for the isotropic bending case, longer remains very close to the reference continuum solution.
To further support the conclusion, that the introduced structural plate model shows very satisfying results, we consider a more complicated loading scenario. The loading is applied quasi-statically, the bending moment is now not only increasing monotonously in time, but its direction is rotating, such that the assumption of proportional loading no longer holds and the continuum solution is not directly available. The tensor of bending moment changes now in time according to the law

M(t) = M(t)(cos 2 (t)e
x e x + cos (t) sin (t)(e x e y + e y e x ) + sin 2 (t)e y e y ), with the magnitude of the bending moment M and its rotation angle varying linear in time: the dimensionless time variable changes from 0 to 1. Plasticity starts with the bending moment M e * in x direction, and in the end of the simulation we reach 95% of the maximal possible bending moment M p * in y direction. We begin the analysis with the structural model. Equilibrating the externally imposed moments (62) with the elastic response (4) and furthermore considering the additive decomposition of the total curvature (6), we obtain the components of the tensor of total curvatures: We substitute these relations in the governing differential Equations (11), (12), and (14). The differential equations for the plastic curvature components and the dissipative work follow after elimination of the plastic consistency parameter : We restrict the consideration to the hardening function (43), identified for the case of uni-axial bending, and numerically integrate the obtained system of differential equations. Using (64), we find the time histories of the components of the total curvature tensor, which are plotted in Figure 3 in comparison to the purely uni-axial case = ∕2, that is, when the bending moments acts along the y-axis all the time: M = M(t)e y e y .
The final values at t = 1 with = ∕2, which correspond to the same bending moment M = 0.95M p * e y e y , are clearly different. From these results we see, that the time history of the applied bending moment affects the current elastic-plastic F I G U R E 3 Comparison of time histories of curvature components for rotating bending moment versus uni-axial moment response of the plate. A few other observations can be made. In the very beginning at t = 0 both solutions are yet purely elastic and rotated with respect to each other, in the same way as the applied moments. Here we also notice that the ratio between both diagonal components of the elastic curvature is exactly . The non-diagonal curvature component xy becomes nonzero only when the moment is rotating.
Ultimately, we seek to validate the results against the continuum model for the case of rotating moment. In the absence of an analytical solution, we need to replace the integral for the total bending moment (19) by a sum using a quadrature formula with discrete integration points. Each of the points contributes to M, and the local strain components in the plane of the plate are determined by the total curvature similar to (17). This makes the moment-driven problem formulation complicated, because the elastic-plastic deformation processes in all the integration points are coupled. We simplify our task and use the previously obtained time history (t) from the solution with the structural theory as an input for the continuum model, which makes the computation straightforward. The plate kinematics provides us with the plane part of the three-dimensional strain tensor − z. Presuming the plane stress state in the entire through-the-thickness element and using the yield condition in the form (48), we integrate the local stress-strain relation over time in each integration point, find stresses (z, t) and then compute the time history of the resulting moments M(t) by integration of (19) using five equidistant integration points in the upper half of the through-the-thickness element 0 ≤ z ≤ h∕2. Time histories of the components of the computed moment tensor are compared against the originally imposed one (62) in Figure 4. The observed correspondence convinces us, that the developed structural theory consistently describes the constitutive behavior of the significantly more complicated continuum plate model in the broad range of loading conditions.

STRUCTURAL PLASTICITY MODEL IN A FINITE ELEMENT SIMULATION
We approximate the deflection field w(x, y) using C 1 -continuous four-node Bogner-Fox-Schmitt finite elements with bi-cubic Hermitian shape functions. 24,25 The functional of the elastic strain energy of the isotropic Kirchhoff plate which comprises the contribution of each finite element with the area Ω el , becomes thus a function of the entire vector of nodal variables W. In each node i we have four degrees of freedom, namely the deflection w i , its first derivatives ( 1 w) i , ( 2 w) i as well as the mixed second-order derivatives ( 1 2 w) i , in which F I G U R E 4 Comparison of time histories of moment components of continuum versus structural theory stands for the derivative with respect to the local coordinate on the element q , = 1, 2. Rectangular elements are used in the following, such that the Cartesian coordinate x depends only on the local coordinate q 1 , and y only on q 2 . The variable element size in the below considered examples and the requirement regarding continuity of the deflection gradient ∇w imply, however, that the function x(q 1 ) is nonlinear, such that the geometry is also smoothly approximated across the element boundaries. We seek the state of equilibrium by minimizing the total energy of the system with the potential energy of external distributed forces the load factor grows incrementally from 0 to 1. In the considered geometrically linear setting the energy U Σ is quadratic in the kinematic variables W. For a given distribution of plastic curvatures p (x, y), the energy is minimized by solving a system of linear equations for W: Here K is the stiffness matrix of the elastic plate, the force vector F e is determined by the externally applied loading, and the vector F p depends on the distribution of the plastic curvature, because of e = − p in (66). The integrals (67) are evaluated numerically using 3 × 3 Gaussian quadrature rule. Plastic curvature values need to be available in each integration point of the model.
As soon as the yield criterion is violated in at least one of the integration points, the rate Equations (11)-(14) need to be integrated over time. For sufficiently small increments of the load factor , the yield criterion function of the new moment tensor f (M) exceeds the current value of k 2 M by just a small value, and we shall use linearized approximation of both terms in order to accurately integrate the evolution law over time. The update of the local plastic variables p and A p within a load increment is performed such, that both sides of the consistency condition (71) become equal for the new values. The flow rule (11) closes the problem, making the solution for the local increments Δ p and ΔA p unique. This leads us to the following formulation of the classical elastic predictor -plastic corrector return mapping algorithm 26,27 for the structural plasticity model.
From the differential relation for the dissipative work (12) we find Demanding the consistency condition f + Δf = k 2 M + Δk 2 M to be fulfilled after the update, we write Here the value of f is to be evaluated for the new deformed state with updated kinematic degrees of freedom W. Its derivative with respect to the consistency parameter follows by first substituting the flow rule (11) into the additive decomposition of the curvature tensor (6) considering the total curvature = const, then finding the time derivative of the moment tensorṀ from the elastic law (4), and finally substituting the result into (9) and dividing bẏ. In Cartesian components the result of this computation reads Having computed Δ , we immediately obtain the increment of the dissipative work from (72) and that of the plastic curvature from Clearly, the update of the plastic variables changes the internal moments in the plate and results in the violation of the equilibrium equation. Updating the kinematic degrees of freedom and thus the deformations according to (70) once more, we again need to adjust the plastic variables, which results into a kind of fixed-point iterations, see Reference 28. The iterations will in most cases converge to a state, in which all conditions are fulfilled and which is thus the solution of the problem for the current value of the load factor . The next load increment shall be processed.
The algorithm is simple and each iteration is efficient, as the constant elastic stiffness matrix K can be decomposed in advance. This advantage is, however, outweighed by the slow (linear) convergence rate of the iterations. The asymptotic quadratic rate of convergence of the Newton method is achieved by using the elastoplastic tangent stiffness matrix K t , which is no longer a matrix of the quadratic form of the strain energy, but rather originates from the equilibrium condition in form of the principle of virtual work. Thus, the virtual work of external forces is a linear form the virtual work of internal forces reads and the principle of virtual work A e + A i = 0 leads to the equilibrium equation The variations of the nodal variables W determine the virtual deflections w and variations of total curvatures via the finite element kinematics. The coefficients of the linear form of the virtual work of internal forces are Here the derivative of the curvature tensor with respect to the degree of freedom W j follows directly from the kinematics and does not depend on W in the geometrically linear model. In the absence of active plastic flow we obtain F i = F p − KW and further again (70). The Newton method for the equilibrium Equation (78) features increments of nodal variables ΔW, which result from the linear system of equations with the system matrix being the external forces are constant. Using (79), we compute the components of the tangent stiffness matrix: Finding the derivative of the tensor of moments is easy in the elastic zones, because here the rates of the total curvature and of the elastic curvature are equal. The symmetric fourth rank tensor of the derivative of the moment with respect to e follows directly from the elasticity law (4). In the points with active plastic flow, the derivative of the moment tensor in (82) must take into account the nonvanishing rate of plastic curvature, such that d e = d − d p , and Similar to the incremental form (73) and (75) we write Note that the derivative of the yield criterion function is computed for constant plastic curvature p . As it is always the case when the associated flow rule is used, the tangent stiffness matrix is symmetric: K t ij = K t ji , which greatly simplifies the implementation. This symmetry follows directly from the definition (82) and formulas (84), (85). In the numerical experiments below, the Newton iterations with the updated matrix K t converge rapidly as soon as the set of integration points with active plastic flow is no longer changing significantly between the iterations.
The discussed backward time integration scheme produces accurate results for small load increments. In the numerical examples below we applied an adaptive refinement scheme. The load increments Δ are reduced, if in the course of Newton iterations the intermediate value of f in (73) exceeded the current level of k 2 M by more than 50%, which guarantees accurate solution. From the point of view of implementation, it is efficient to store in each integration point the current values of and k M along with the primary plastic variables and to keep them constant during the Newton iterations for a given load increment. This greatly reduces the number of calls to the computationally expensive numerical solver for the Equation (59), which is relevant if material hardening is taken into account.
For the sake of comparison, in the following we also present results of computations with the continuum plasticity model, which features five equidistant integration points in the upper half of the plate 0 ≤ z ≤ h∕2. In each integration point the plane stress plasticity equations are integrated using the conventional return mapping algorithm with the tangent stiffness matrix according to Reference 26.

BENCHMARK RESULTS FOR A SQUARE PLATE
We begin by considering the conventional numerical example of a simply supported square plate under uniformly distributed loading q = const, without material hardening (see Reference 10). Upon incrementally increasing the load, we observe that the plastic deformation first appears in the corners of the plate. Because the bending moment tensor near the edges of the plate cannot be isotropic we expect the hardening function k M , which is identified by uni-axial bending experiment, to be more appropriate in this case. Motivated by this consideration, we used the expression (44) for the below simulations with the structural plate model. The model parameters are according to Reference 10 and are provided in Table 2; = b is the length and the width of the plate. Prior to comparing the results of the general structural plate model to the other simulation, we conduct a convergence study by virtue of computing the inelastic response in form of the displacement of the center point of the plate w mid for growing load values q. The results are provided in Figure 5 for various meshes.
Each marker in this load-displacement diagram corresponds to a certain load level q; note the uneven distribution of load increments, which results from the limitation of the maximal plastic flow for the sake of accuracy of the numerical integration. We observe that the deflections for the same load levels become barely distinguishable for the two finer finite element discretization levels even in the higher range of q when the structure approaches the plastic collapse. Therefore, We proceed by comparing the obtained load-displacement diagram to a full continuum solution with a solid model of the three-dimensional body of the plate computed in ABAQUS and to a continuum plasticity model of the plate, which features five equidistant integration points in the upper half of the plate 0 ≤ z ≤ h∕2 and is handled by the classical plasticity formulation with the return mapping algorithm according to Reference 26. Since the latter model is handled within the same finite element discretization which is used for the structural plate model, the mesh and the type of the elements are identical. As for the full continuum ABAQUS model, a uniform mesh of 100 × 100 C3D20 (20-noded quadratic brick) elements in the plane and four elements over the thickness is used.
In Figure 6, very good compliance between the new general structural plate model and classical continuum approaches is seen. Comparing to the solution results available in Reference 10 (dots in Figure 6), we notice certain differences, although both solutions are obtained without the consideration of material hardening. Despite the use of the shear-deformable Reissner-Mindlin plate theory by the authors of Reference 10 and the relatively high thickness of the plate, the successful comparison against the full three-dimensional solution convinces us that the new structural plate model is indeed sufficiently accurate. The observed lower accuracy of the results of Reference 10 is explained by the fact, that this model does not take into account the spreading of plasticity through the thickness. Nevertheless, it is capable of predicting the limit loads. The lower and upper bounds for the estimation of the limit load parameter (plastic collapse) are analytically computed in Reference 10 according to the results of a limit analysis by the application of the yield-line theory, see Lubliner: 8 Here = b for the square plate (see Table 2) and M p * follows from (26). From Figure 6 we conclude that like in Reference 10 our limit load (around q = 0.61) also falls within the given analytically estimated interval as it should.
In Figure 7 we demonstrate the results of a simple convergence study regarding the effect of the number of integration points over the thickness in the continuum model. It justifies the necessity to use five equidistant integration points over the upper half of the plate. The computationally cheaper models with just 2 or 3 integration points result into a considerably less-accurate load-displacement diagram, which underlines the relevance of the present study.
The actual gain in the efficiency of the structural model is demonstrated by the comparison of computation times, needed to obtain the results of Figure 6. From Table 3 we see, that the stress resultant model is almost twice as efficient as the traditional continuum one. Similar relation of the efficiency was observed for further benchmark examples as well.
We proceed with the comparison of the growth of the plastic zones between the structural plate model and the continuum plate model.
In Figure 8  (C) q = 0.55 lines for structural plate) enclose areas with nontrivial plastic strains in the integration points. We observe that a fifth plastic zone appears in the center of the plate as q is growing in addition to the four initiated in the corners. At some load levels they all merge into a single large plastic zone. Again, good agreement between the structural mechanics and the continuum plasticity models is observed. We conclude this section by analyzing a loading-reverse loading-loading scenario for the here considered simply supported square plate under uniformly distributed loading. Although it is theoretically clear, that the presented structural model cannot accurately reflect reverse plasticity because of the purely isotropic hardening, the experimental demonstration sheds more light onto the issue. In Figure 9 we compare load-displacement diagrams with nonmonotonous loading history for the continuum and the structural mechanics models.
Both solutions are very close during the initial loading. The maximal force levels in both models are slightly different and chosen such, that the middle point deflections in the end of the initial loading stage are the same and equal w mid = 5 ⋅ 10 −4 . Such displacement-controlled loading results in a more informative diagram.
During the reverse loading we steadily change the direction of the distributed force, while its magnitude in the end of this stage is the same. The observed differences are explained by the Bauschinger effect in the continuum plate model, in which yielding occurs earlier than in the structural plate model with just the isotropic hardening. The symmetric hysteresis loop of the continuum model does better correspond to the theoretical expectations. 8,29 Nevertheless, results presented in Figure 4 support the conclusion, that the range of applicability of the structural model reaches far beyond the purely proportional and monotonous loading scenarios.

BENCHMARK RESULTS FOR A RECTANGULAR PLATE
Another appropriate benchmark problem with simulation results available in the open literature is a rectangular plate under uniformly distributed loading with ideal elastic-plastic material, considered by Dujc and Brank 12 with the account of gradual spreading of plasticity through the thickness. The model used in Reference 12 features an adapted Ilyushin yield criterion with the time-varying yield moment m * 0 ( ) = (t)m 0 as proposed earlier by Crisfield. 11 The time evolution law for the field variable (t) is designed such, that it matches the uni-axial bending experiment. This allows capturing the mentioned effect of the development of plastic zones in a through-the-thickness element and results in a good correspondence with continuum mechanics solutions. The degree of this agreement depends, however, on the boundary conditions, as their choice determines the locations and shapes of the plastic zones in the plate. Therefore we chose the case with simply supported edges, as it results into a slightly worse agreement between the stress resultant and ABAQUS solutions in Reference 12. The parameters of the benchmark example are summarized in Table 4.
We discretized the domain with 60 × 40 finite elements, which formally corresponds to the discretization used in Reference 12-although with different shape functions. In Figure 10 we compare the results of our structural mechanics model to the one from Reference 12 (obtained by sampling the graphical data of the published paper). We also plot the F I G U R E 9 Uniform loading-displacement for the loading-reverse loading-loading scenario: structural plate versus continuum plate solution of the plate model with 60 × 40 S4R (quadrilateral shell) elements with 10 integration points over thickness, computed in ABAQUS for the sake of reference.
We see, that the stress resultant plate model of Reference 12 tends to overestimate the deformation in the earlier plastic phase and to underestimate it in the latter plastic state. The slight nonsmoothness of the respective curve is due to the time integration strategy for the parameter (t), which is considered constant within a time step. While the model 12 shows much better correspondence to the reference ABAQUS solution for a clamped plate, the accuracy of the present approach appears to be essentially less affected by the particular choice of the boundary conditions and loading.

BENCHMARK RESULTS FOR AN ELONGATED PLATE WITH SELF-EQUILIBRATED LINE LOADING
In the last example problem we address an elongated rectangular plate ( > b) with loading, which resembles the force distribution occurring during the roll forming process. The left boundary of the cantilever plate is clamped and all other boundaries are free. The self-equilibrated loading is distributed along the line in the width direction in the middle of the plate: q(x, y) = P(y) (x − ∕2) with being the Dirac impulse and the load intensity per unit length of the line is p is the amplitude of the load. Simulation parameters are provided in Table 5 and the distribution of the load along the width of the plate is depicted in Figure 11.
Due to the fact that the considered plate problem shows similarity to a cantilever beam, again the hardening function in the form (44), which is identified for uni-axial bending, is used for the structural plate model.
In the simulations we consider a mesh with 20 elements over the width, while the element size in the length direction is variable, such that finer discretization is reached in the middle of the plate, where the load is acting, see Figure   Parameter b  15. We increased the loading factor in the model, until the specified load amplitude from Table 5 is reached. Afterwards we fully unload the plate, in order to evaluate its new shape as a result of the residual plastic strains. The corresponding deformed shapes of the structural plate model and of the classical continuum plasticity model are shown in Figure 12.
At the edge x = 0 the plate is clamped and no deformation is visible there. The plotted configuration with the larger deflection downwards at the right end corresponds to the residual deformation of the structural plate. This slight difference in the residual deflections far away from the location of the loading should be regarded as insignificant owing to the two reasons: • it is the shape of the residual deformation of the plate, at the location of action of the rolling forces, which is of actual practical importance for modeling the rolling process and • minimal differences in the deformed state in the middle of the plate, where plastic deformations are essential, would result into noticeable change of the deflection of its right end because of the long lever arm.
The comparison of both solutions is better visible from the distributions of the deflection in the length direction w(x, 0) provided in Figure 13 and in the width direction in the middle of the plate w( ∕2, y) plotted in Figure 14, which further underlines the similarity between the computed deformed shapes within the middle cross-section as a result of the rolling action.
Furthermore, we again highlight the development of the plastic zones, whose configurations for three levels of the load amplitude p are depicted in Figure 15.
Three plastic zones are observable at lower values of p, while they later merge into a single one with two small elastic "islands" at the middle transverse line near the points, where the line load intensity P(y) vanishes.
We conclude the analysis by discussing the results of the structural plasticity model with the account for material hardening, which is important to consider in metal forming processes. Making use of the Swift law (56) with the strain hardening exponent value n = 0.26, we obtain deflections, plotted by dotted lines in Figures 13 and 14. While the solution F I G U R E 11 Self-equilibrated line load, which acts along the width, at the center of the plate F I G U R E 12 Residual deformation of the structural plasticity and of the continuum plasticity plate models; note the scaling of the deflections by the factor of 2 ⋅ 10 3 , for the sake of better visualization F I G U R E 13 Distribution of the residual transverse deflection w(x, 0) along the center line of the plate in the length direction, computed for three models of plastic behavior F I G U R E 14 Distribution of the residual transverse deflection w( ∕2, y) along the middle transverse line of the plate in the width direction, computed for three models of plastic behavior does not change qualitatively, the magnitude of the residual deformations expectedly becomes smaller. The shapes of the plastic zones for the same load values change slightly.

CONCLUSIONS
Presently, simulations of elastic-plastic response of structural members mostly make use of three-dimensional continuum plasticity laws. Already published attempts to handle plasticity solely on the structural mechanics level and thus to obtain a formulation based on stress resultants only, did not really make it into computational practice because of the issues with reliability, convenience of use, and uncertainties regarding the parameters of the model. In this paper we introduce a novel approach, which we expect to be convenient and efficient for practical applications because of its simplicity, consistency and sound theoretical basis. Thereby, the treatment of the classical inelastic constitutive relations of continuum mechanics in several integration points along the thickness becomes obsolete. This has a remarkable positive impact on the computational efficiency of the simulation. Thus, the computations for above benchmark tests took approximately half the CPU time with the new structural mechanics formulation. We expect the proposed model to be particularly efficient in the field of roll forming simulations, which in their own turn need to be fast for the practically relevant purposes of efficient development of the production process as well as of model-based controller design. The basic idea of the proposed structural mechanics plasticity model rests on the identification of a hardening function in the equations featuring curvatures and bending moments. We perform this identification using available analytic solutions for simple loading cases, such as uni-axial or isotropic bending. We restrict the consideration to the isotropic hardening law, because reverse plasticity is not common in roll forming processes. The formulation is presented for both the ideal elastic-plastic material behavior, as well as for isotropic material hardening according to a power law, which is conventionally used in the field of material forming. The finite element structural mechanics model is implemented in the classical form of a return mapping algorithm with tangent stiffness matrix. The proposed approach is validated by comparison against both the full three-dimensional solutions of a benchmark problem as well as against solutions, obtained with the plate theory using integration over the thickness. The comparison results show good accuracy for the considered class of benchmark examples. For future research, we aim at further extending the formulation towards a geometrically nonlinear setting with nonmaterial Eulerian-Lagrangian kinematic description in the spirit of References 4,5. Along with the account for contact conditions, this would allow for realistic simulations of the roll forming process.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.