A Constitutive Model for Stud Connection in Composite Structures

: The complexity of finite element analysis for composite structures can be significantly reduced by representing the connector and adjacent concrete as a macroscopic element. Nevertheless, the prevailing macroscopic models for shear connections predominantly employ nonlinear elastic theory. This approach introduces inaccuracies in estimating structural stiffness and load-bearing capabilities, primarily due to its inability to precisely capture the cumulative effects of plastic damage. In response, this study introduces a novel macroscopic elastoplastic model grounded in plasticity theory, aimed at accurately characterizing the nonlinear behavior of stud connections subjected to concurrent shear and tensile forces. This paper meticulously delineates the implementation of the elastoplastic constitutive model using the backward Euler method for numerical integration. It further articulates the derivation of the consistent tangent stiffness, which aligns with the convergence efficiency of the Newton–Raphson iterative approach. The computation of the element stiffness matrix for a two-node element is executed via the governing equation inherent to the finite element method. An exemplar macroelement test conducted in ABAQUS affirms the implicit backward Euler scheme’s stability and consistency across varying tolerances. Validation of the elastoplastic model against empirical test outcomes corroborates its efficacy, demonstrating the model’s precision in predicting the load–displacement behavior of stud connections under the influence of shear and tensile forces


Introduction
In the domain of composite structural design, accurate forecasts of shear connections' performance are critically important.These connections serve as a means for transferring both axial and shear loads between distinct components of composite structures, and should they fail, the resulting consequences may prove to be catastrophic [1][2][3][4].When a composite beam is only subject to bending moments and shearing forces, the inherent force within the shear connections is predominantly shearing [5][6][7][8][9].Accordingly, past research has tended to concentrate on the shear strength of these connections under such circumstances [10][11][12][13][14].However, when subjected to torque or localized loadings, these shear connections experience an amalgam of shear and tensile forces [15][16][17][18].Evidence of the coupling effect within stud connections concerning tensile and shear strength was presented in investigations by Lin et al. [19].Studies by Tan et al. [20] examined the capabilities of demountable steel-concrete connectors under the simultaneous presence of shear and tensile loads.This study highlighted a noteworthy reduction in shear resistance when tensile forces were introduced, signifying the complex interplay between the two types of forces.Practical applications of stud shear connections in both homogenous concrete slabs and composite slabs were explored by Shen and Chung [21].Their findings revealed substantial variations in ductility at minimal and substantial deformation levels.Collectively, these investigations stress the imperative nature of factoring both shear and tensile forces into stud connection designs.By doing so, the connections are optimized for efficient load transmissions, thereby decreasing potential safety hazards.Thus, an investigation into the behavior of shear connections under the concurrent influence of shear and tensile loading certainly holds great potential for enhancing safety provisions for composite structures.
Throughout regular usage, failures of shear connections occur comparatively infrequently.Consequently, the impact of shear connections on the load-bearing capacity of composite structures mainly manifests through the load-slip relationship [19,[22][23][24].Prior research attests to the nonlinear elastoplastic behavior that shear connections display during the loading process [25,26].In prior analyses, a nonlinear elastic constitutive model was employed, primarily due to the absence of an elastoplastic constitutive model [27,28].Nonetheless, the utilization of a nonlinear elastic constitutive model can introduce inaccuracies in the examination of structural rigidity and load-bearing capability.It fails in accurately portraying the internal force redistribution and cumulative accrual of plastic damage under realistic circumstances [29,30].Therefore, the development of the constitutive elastoplastic model of shear connections under the simultaneous influence of shear and tensile loads becomes of paramount importance.
Plasticity theory has long found its application in the constitutive modeling of materials such as metals and soils [31][32][33].Current trends see the theory being increasingly leveraged to stimulate the multiaxial loading behavior of shallow foundations or oil pipelines on a macroscopic level [34][35][36][37][38].These constitutive models offer a holistic representation of foundation behavior, establishing a direct linkage between the load and the corresponding displacement.Since the foundation response is articulated in consistent terms as structural mechanics, the load-displacement relationship can be seamlessly integrated into the numerical structural analysis.In a similar vein, assorted phenomenological models are routinely employed to recreate the macroscopic nonlinear behavior of complex materials and structures [39,40].These macro-level constitutive models bypass any special interface or contact, providing a practical approach for the interactive modeling of foundations and the overall structures.
Borrowing from the macroscopic constitutive models of shallow foundations and oil pipelines, this paper proffers an exploration of the elastoplastic constitutive model of stud connections facing combined shear and tensile loads.Through experimental results and finite element analysis, coupled with an extensive literature review, a formula reflecting the load-plastic displacement interaction of stud connections under taxing load conditions is proposed.A formulation for the final failure surface under the simultaneous effects of shear and tensile loads is set down, which births the surface expression for an elastoplastic model.Based on experimental results and finite element analysis, alongside an extensive review of the literature, an expression for the load-plastic displacement of stud connections under complex loading conditions is proposed.An expression for the ultimate failure surface under the simultaneous impacts of shear and tensile loads is formulated, yielding the surface expression for an elastoplastic model.An isotropic hardening rule specific to stud connections is proposed, followed by the introduction of a flow rule and a plastic potential to accommodate the plastic deformation characteristics.By utilizing an implicit integration algorithm grounded on a backward Euler method, updating the force vector and capturing plastic deformation behavior becomes achievable in the overall constitutive model.This facilitates the realization of an elastoplastic model for stud connections and prompts its numerical implementation.By developing user-defined subroutines, the elastoplastic model is applied to the calculation of push-out tests to verify the performance of the constitutive model of stud connections.

Outline of the Model
Positioned at the intersection between the steel beam and concrete slab, the stud connector is subjected to diverse permutations of vertical and horizontal movements under the simultaneous impacts of shear and tension, subsequently engendering a complex stress state.The meticulously defined macroscopic elastoplastic constitutive model successfully captures the response of the stud connection through defining the relationship between the force vector, denoted as F, and the relative interface displacement vector, denoted as U.
This relationship, designated in-plane, is hereinafter referred to as the load-displacement correlation.The definitions of F and U imply that the composite response of both the stud connector and the adjoining reinforced concrete are implicitly encapsulated within the model.This observation further dictates that any variations-be they material properties or geometrical nuances-within either constituent of the stud connection exert a significant influence on the parameters of this model.The directions of the net positive force and the consequent relative displacement at the steel-concrete interface resulting from the forward force can be inspected in Figure 1.Given this backdrop, it is indeed pertinent to define F and U in the following manner: Attention should be directed to the fact that within this paper, variables presented in bold typeface consistently signify vectors or matrices.In the context of capturing the nonlinear behavior of stud connections under the synergistic influence of shear and tensile loads, the macroscopic elastoplastic model is primarily composed of the following elements: (1) An empirical formulation of the yield surface within the force plane, labeled as (N, V).This is achieved by delineating the contour of the yield surface, which is informed by the shape of the failure envelope; (2) A strain hardening scheme that elucidates the functional interrelation between the magnitude of the yield surface and the plastic deformation.The expansion of the yield surface is a function of the evolution of plastic displacement U p (where U p = u p v p T ).
The extent to which the yield surface expands is quantified by the shear force V 0 that corresponds to a given plastic displacement; (3) An aptly chosen flow rule that enables the prediction of increments in plastic displacement, denoted as dU p , which results from the yield surface's dilation or constriction.(4) An elastic load-displacement relationship is employed to delineate the load-displacement behavior that occurs within the confines of the yield surface.
The acquisition of the components integral to the above-outlined elastoplastic model mandates the undertaking of an adaptation process.This process involves fitting expressions, representative of these components, to data derived either experimentally or through the deployment of finite element simulations.This article espouses the utilization of three distinct loading paths, which serve to offer ancillary data support pivotal to the fabrication of the elastoplastic model.These chosen trajectories are illustrated in Figure 2: (1) Load path 1: shear direction loading under different drawing force levels.The shape enclosed by the extremities of the load path may be posited as depicting the trajectory of the yield surface which is apt for the elastoplastic model.This load path is instrumental in deducing both the hardening rule and the flow principle inherent to the elastoplastic material model; (2) Load path 2: radial loading.These loading trajectories inform the hardening principle as well as the flow rule characteristic of the elastoplastic model.The instance of pure shear loading does indeed represent a specific case within this load path, from which the shear-slip hardening principle applicable to the stud connection can be extrapolated; (3) Load path 3: elastic stiffness test.By introducing minimal displacements in both planar directions while concurrently unloading the stud connection, it is feasible to establish an estimation of the elastic stiffness coefficient.
Load path 2

V N
Load path 1 Load path 3 To this end, a series of push-out tests on stud connections were conducted, using the test setup depicted in Figure 3.The employed configuration incorporated a hydraulic jack positioned between the concrete slabs designated to impose tensile stress upon the stud connections, extending along with a separate jack linked to the reaction frame to exert shear forces.The mutable parameters within the push-out test specimens included the elastic modulus of the concrete, the diameter of the stud, and the intensity of the applied tensile force onto the stud.The shear-slip curve under different tensile forces was thus obtained through these tests for the stud connections.Subsequently, a finite element modeling approach was implemented to simulate the push-out tests under shear load for the stud connections, using the finite element model shown in Figure 4.The finite element method empowered us to imitate a broader spectrum of loading paths and engage in parametric analysis.This occasioned an intensive exploration into the ramifications of variables such as the length and diameter of the stud, the elastic modulus of both the concrete slab and the stud, and the volume ratio of reinforcement.Whilst a summary of the elastoplastic model is present, a more comprehensive breakdown of the push-out tests and the finite element models is available in the cited work [41].

Elastic Behavior
The elastic stiffness inherent in the stud connection constitutes a critical facet of its associated elastoplastic model.In the pursuit of computing any load or displacement increment occurring within the yield surface, it is imperative to furnish a characterization of the stud connection's elastic response.The elastic relationship between the increment of force ∆F and the corresponding increment of elastic displacements ∆U e can be expressed as follows: where K e represents the elastic stiffness matrix, k 11 and k 22 denote the elastic tensile and shear stiffness, respectively, and k 12 = k 21 is the coupling stiffness.The elastic stiffness matrix K e is obtained by fitting the result of the finite element simulation of the stud connections under load path 3 [41], and it can be expressed as follows: where d s , E c , and E s are the stud diameter, the elastic modulus of the concrete, and the elastic modulus of the stud, respectively.

Yield Surface
The demarcation between the purely elastic and the combined elastic-plastic states is referred to as the yield surface.Specifically identifying the first transition from an elastic to an elastoplastic state, the defined yield surface can be called the initial yield surface for the stud connection.As the loading process progresses, this initial yield surface undergoes a constant expansion, generating what can then be described as the subsequent yield surface.Bringing this process to conclusion, when the connection attains its point of failure, the yield surface present at that juncture is coined the failure surface.Figure 5 provides a graphic representation of these stages, displaying schematic outlines of the initial yield surface, the subsequent one, and, ultimately, the failure surface.Regarding the shape of the failure surface, established experimental data and data from finite element simulations [42][43][44][45] have collectively affirmed that it approximates to the envelope of a convex arch.On this evidence, we can represent the failure surface through an equation that captures the form of such a convex envelope.
where N 0max and V 0max are the radius of the horizontal and vertical axes of the failure surface, respectively, and β 1 is the aspect ratio defining the surface shape.The initial path of the yield surface in the (N, V) load plane can be obtained by shrinking the shape of the failure surface as follows: In order to derive the mathematical expressions for both the failure surface and the initial yield surface, as denoted in Equations ( 5) and ( 6), it is imperative to ascertain the formulae for the characteristic parameters N 0max , V 0max , V 0i , and β 1 .The values of these parameters are deduced by fitting them to the empirical outcomes obtained from push-out tests, as well as the results from finite element simulations corresponding to load path 1, as documented in the study by Qin [41]: Presuming that the morphology of the subsequent yield surface mirrors that of both the initial yield surface and the failure surface, we can formulate the subsequent yield surface.This is achieved by implementing an expansion of the initial yield surface, specified as follows: The size of the yield surface is represented by V 0 , the shape of the yield surface is determined by n 0 (n 0 = N 0max /V 0max ), and the centroid of the yield surface is (0, 0).

Hardening Law
Building on the hypothesis that an augmentation in plastic displacement within the shear dimension contributes to the expansion of the yield surface, the size of the yield surface when the system is in a pure shear state is dictated by the principles of isotropic hardening.The isotropic hardening law delineates the interdependence of the shear force V 0 and the plastic slip v p pertinent to the shear direction in the stud connectors.To encapsulate this, we draw on the shear-slip correlation derived from the insights of push-out tests complemented by finite element analysis simulations.The equation representing this correlation is: Further, the relation between the shear force V 0 and the plastic slip v p in the shear direction of the stud connection during the elastic-plastic stage is where v p represents the cumulative outcome of integrating the incremental plastic displacements dv p along the shear vector.Drawing from the empirically derived shear-slip response curves of the push-out tests under pure shear conditions, as documented by Qin [41], a functional relationship is established between the shear force V 0 and the total slip v.To isolate the plastic component, one subtracts the elastic part from the total slip v, thereby obtaining the plastic slip v p .By differentiating the shear force V 0 with regard to the plastic slip v p , we arrive at a differential equation that articulates the dynamic interplay between the shear force V 0 and the plastic slip v p within the realm of the shear direction: Although Equation ( 11) might be sufficient to fit the shear-slip curve, it does not account for complex loading conditions that entail plastic displacement in both the shear and tensile directions simultaneously.To account for the influence of plastic displacement in the tensile direction on the size of the yield surface V 0 , an alternative methodology is necessary.To maintain consistency with the hardening law for both the pure shear and shear-tensile states, a new parameter called the plastic displacement strength, denoted as ∆ p , is introduced to replace v p in the formula when the stud connection experiences plastic displacement in both the shear and tensile directions, namely The plastic displacement strength ∆ p is a function of the plastic displacement v p in the shear direction and the plastic displacement u p in the tensile direction, which can be obtained by integrating the plastic displacement increments dv p and du p , respectively.It is crucial to ensure the hardening law observed in the shear-tensile states aligns seamlessly with that strictly within the pure shear domain.To discern the dimensions of the yield surface, one must compute the value of V 0 as derived from the vertical plastic displacement v p0 during a pure shear state.Under shear-tensile conditions, ∆ p equates to v p0 .This conclusion can be reached by replacing the corresponding values of v p and u p at any point on the yield surface with ∆ p .By adhering to these guiding principles, one can acquire the optimally fitted value for Equation (12) as follows:

Flow Rule
The increase in plastic displacement at yield defines an appropriate flow rule.The flow rule within the two-dimensional (N, V) load plane is informed by load paths 1 and 2, and a non-associated flow criterion espoused by the elastoplastic stud connection model is exhibited [41].To model the interplay between load and incremental displacement, it becomes requisite to charter an expression for the plastic potential g, which retains its identity distinct from the yield surface.Herein, our focus shifts to presenting a formulation for the plastic potential, one that mimics the equation for the yield surface closely, but introduces modifications at the apex and shape of the plastic potential surface.This adaptation is enabled using plastic potential parameters V ′ 0 , α n , and β 2 .Consequently, the equation delineating the plastic potential can be articulated as follows: where V ′ 0 is the maximum shear value of the plastic potential at the present time, which is similar to V 0 in the expression of the yield surface.By regulating the parameters α n and β 2 , one can subtly alter the form of the plastic potential surface.
The ideal values of α n fluctuate in accordance with the varying ratios of plastic displacement iterations dv p /du p .To uphold the flow rule's consistency across an assortment of dv p /du p values, this nuance must be carefully addressed.To accurately emulate this dynamic, an escalation in α n is required as dv p /du p increases, highlighting an intensified degree of non-associated flow.One can chart the fluctuation pattern of α n employing hyperbolic functions.In this context, the plastic displacement iteration dv p , as observed in the shear direction, is rendered normalized by the corresponding plastic displacement iteration du p within the tensile direction: where α n0 represents the parameter value of α n with no increase in plastic displacement, α n∞ represents the limiting value attainable by the parameter, and k ′ n represents the rate at which α n attains its limiting value α n∞ .Figure 6 depicts the evolution of α n .The impact of fluctuating α n values on the theoretical prediction curve is depicted in Figure 7: During the initial timeframe at point A, where α n = α n0 , a compelling proximity is observed between the experimental and finite element results in conjunction with the theoretical prediction curve.With load application, the experimental and finite element outcomes progressively gravitate towards the theoretical prediction curve, under circumstances of α n = α n∞ .The pace of this transition is dictated by the expansion rate of the plastic potential surface, which, in turn, is determined by k ′ n .
The relationship between α n and dv p /du p .
Equation ( 14) articulates the theoretical variations in the plastic displacement within both the shear and tensile directions: where λ can be derived from the consistency condition intrinsic to the strain hardening law.The enhancement observed in the plastic displacement ratio can be represented through an expression relevant to the current state: As depicted in Figure 8, the angle θ vu (theory) (= arctan v p /u p (theory) ) signifies the direction of the theoretical plastic displacement increment, as calculated by Equation (18).On the other hand, θ vu (exp) (= arctan v p /u p (exp) ) characterizes the direction of the plastic displacement increment, as derived from experimental and finite element simulation results.In seeking a fit between θ vu (theory) and θ vu (exp) , the values of the parameters intrinsic to the flow rule emerge as follows: β 2 = 1.69, α n0 = 1.00, α n∞ = 2.40, and k ′ n = 3.00 [41].

7.
The theoretical prediction curves with different values of α n .
To calculate the size of the plastic potential surface (i.e., the size of V ′ 0 ), the parameter x(= V ′ 0 /V 0 ) is defined.The size of x is obtained by the following numerical solution:

Constitutive Integration Algorithm
The elastoplastic behavior of stud connections can be represented through integration with a macroelement, which directly delineates the force-displacement interactions specific to stud connections.This macroelement can be effortlessly integrated into the comprehensive finite element analysis of the structure as a cohesive component.To aid in the customization of macroelement characteristics, Abaqus offers a User Element Subroutine (UEL) [46].The load-displacement curve of the elastoplastic model is encapsulated within this user-defined subroutine.The initial load F n and the plastic state variable V 0n at time t n are passed to the subroutine by the primary analysis program.Upon receiving the displacement increment ∆U n+1 , the subroutine is tasked with computing the subsequent load F n+1 and plastic state variable V 0n+1 at the conclusion of the interval t n+1 using the constitutive integral method, and then relaying these results back to the main program.This symbiotic exchange between the main program and the subroutine is graphically represented in Figure 9. Prior to constitutive integration, it is assumed that the plastic state variables remain unchanged (∆U p n+1 = 0, ∆V 0n+1 = 0) after loading F trial n+1 : This phase is commonly referred to as the elasticity prediction since the trial load F trial n+1 is derived, presuming solely elastic behavior.Within each incremental step, the trial load F trial n+1 is assessed, utilizing either the pure elastic or elastoplastic formulation depending on the load step's termination condition.Five potential scenarios (depicted in Figure 10) are employed to establish the position of the loading state in relation to the currently defined yield surface: Case 1: If f trial < 0, it suggests that the test load resides within the yield surface; hence, an elastic increment should be implemented.The trial load state, F trial n+1 , is deemed the final load state.Thus, F trial n+1 is recognized as the actual force F n+1 , and the value of F trial n+1 should be set as F n+1 : Case 2: If f trial = 0, it signifies that the test load exactly falls on the yield surface.Similarly, an elastic increment is prescribed, and the trial load state F trial n+1 is finalized.Here again, F trial n+1 embodies the actual force F n+1 , and F trial n+1 is assigned to F n+1 .Cases 3, 4, and 5: When f trial > 0, this signifies that the test load is located outside the current yield surface, and, thus, plastic behavior is expected.The constitutive equation will need to integrate the displacement increment: Hence, the primary aim of the constitutive integration is to determine the plastic displacement increment ∆U p n+1 in cases 3, 4, and 5.This research utilizes the backward Euler method for integration [47], which is recognized for being an implicit algorithm.The selection of the backward Euler method is attributed to its computational efficiency and robust theoretical foundation.It is especially beneficial for addressing highly nonlinear issues due to its inherent stability, straightforwardness, and ease of implementation.The implicit algorithm sets out to compute the plastic displacement increment p n+1 and the plastic state variable V 0n+1 , beginning from the trial load's final state F trial n+1 , and adjusts the load state back to the yield surface at t n+1 .Throughout the iteration process, it carries out constitutive integration based on the increment's final state and calculates the load F n+1 at the increment's conclusion.
The plastic multiplier λ and the differential of plastic potential energy ∂g/∂F n+1 in the equation are evaluated using the state at the end of the increment.
As demonstrated in Equation ( 23), an iterative process is needed to navigate the test load back to the yield surface.In the context of the elastoplastic model, F n+1 , U p n+1 , and V 0n+1 should concurrently satisfy the conditions on the updated yield surface.Consequently, alongside Equation ( 23), it becomes crucial to perform integration in accordance with the hardening law to derive the revised plastic state variable: In addition, at the end of the step, the yield condition should be applied: By defining the residual R n+1 , Equations ( 23)-( 25) can be solved by the Newton-Raphson iterative method.In essence, through iteration, the load, displacement, and state variables adhere to the equation of the yield surface, flow rule, hardening law, and consistency condition: The above equation solution of x n+1 = (N n+1 , V n+1 , V 0n+1 , λ n+1 ) can be obtained through iteration, until the residual R n+1 ≤ ϵ tol , and ϵ tol is the allowable error.The iterative procedure is as follows: In the formula, J n+1 is the Jacobian matrix of the above nonlinear Equation, (26), which is a 4 × 4 square matrix.
To improve the clarity and understanding of the integration algorithm process, a detailed flow diagram was constructed, as illustrated in Figure 11.
Process of the integration algorithm.

Continuum and Consistent Tangent Stiffness
The fundamental approach for addressing nonlinear finite element problems is through Newton-Raphson iteration.A crucial step in implementing this iterative method for finite elements is calculating the tangent stiffness matrix, which characterizes the relationship between force and displacement gradients based on the current force and displacement states.

Continuum Tangent Stiffness
The conventional tangent stiffness method usually adopts the continuum elastic-plastic tangent modulus.For the incremental load-displacement relationship of pure elasticity, its expression is The expression of the elastoplastic incremental load-displacement relationship is where K ep is the elastoplastic stiffness matrix, also known as the continuum elastoplastic tangent modulus.
If the load increment remains entirely within the yield surface, the elastic displacement can be determined using Equation (3).Upon yielding, the increment vector of plastic displacement, dU p , is oriented perpendicularly to the plastic potential surface.The magnitude of the plastic displacement increment, dU p , is defined by the scaling factor λ, Whenever there is an increment in load or displacement that crosses the yield surface, the variation in the value of f for the yield function must be zero (i.e., d f = 0) to ensure that the load point, characterized by the normal force N and shear force V, remains on the yield surface.In accordance with the strain hardening rule, a solution for λ must exist signifying that a certain magnitude of plastic displacement is achievable which fulfills the yield criterion where V 0 and dV 0 depend on the plastic displacement strength ∆ p .According to the chain rule and the plastic potential relationship, Equation ( 31) can be written as Therefore, the solution of λ can be obtained from the consistency condition, The continuum elastoplastic tangent modulus K ep can be derived from Equations ( 28)-( 33),

Consistent Tangent Stiffness
The Newton-Raphson iterative method is characterized by a second-order rate of convergence.However, when applied to the computations of the continuum elastic-plastic tangent modulus K ep , the backwards Euler method generally achieves only first-order accuracy.This results in a discrepancy between K ep and the updating format of the pathindependent variables.Consequently, the Newton-Raphson iterative technique cannot guarantee a quadratic convergence rate.To ensure the second-order convergence rate of the tangent stiffness method, it is necessary to devise an updating scheme for path-independent variables that yields first-order precision, aligning it with a consistent elastoplastic tangent modulus D ep [48].Moreover, in finite element analysis, it is crucial to define D ep to derive the Jacobian matrix essential for equilibrium iteration.D ep is defined as where ∆F is the force increment and ∆U is the displacement increment.We know that so Equation ( 35) can be written as where in which I denotes a (2 × 2) identity matrix and 0 denotes a (2 × 2) matrix with zeros.To obtain ∂∆F/∂x n+1 , we reformulate Equation ( 26) as Differentiating Equation (39) yields where the relation d∆F = dF n+1 is used.Therefore, it can be deduced that Inserting Equations ( 38) and (41) into Equation (37) finally yields where A= I 0 is used to obtain a concise form.

Finite Element Formulation
While defining the macroelement through the User Element Subroutine (UEL), the structure of a two-node element is adopted.These two nodes might either coincide or possess a certain distance between them.For the purpose of this explanation, we designate these nodes as i and j.As a result, the displacement increment function vector for the element, denoted as {∆δ} e , can be construed as follows: and the vector of force increment function {∆F} e for the element is The strain-nodal displacement parameter relation of the element is in which L is the length of the element; B is the strain-displacement transformation matrix given by According to the constitutive relations, the stress vector of the element could be written as where A denotes the cross-section area of the element.
Considering the element contribution to the equilibrium of the system, the nodal forces {∆F} e that the element provides are in which V is the volumn of the element, N is the shape function of the element, and {∆f V } is the volume force increment.No volume force is considered in this element, so the equation could be written as One integration point is adopted and the weight factor is one.Substituting Equations ( 45) and (47) into Equation ( 49) results in Thus, the element stiffness is

Application of Integration Algorithm
The efficacy of the elastoplastic constitutive model is substantiated through a numerical illustration.This instance makes use of a singular macroelement force-resultant model (as demonstrated in Figure 12).This study sought to scrutinize the influence of allowable errors, denoted as ϵ tol , amidst the application of the backward Euler integration scheme, particularly with respect to computational precision and efficiency.In this example involving a single macroelement, the stud's diameter was established at 16 mm, with the concrete's compressive strength for a cubic sample measured at 54.06 MPa.The elastic modulus for concrete was set at 35,188 MPa, and for steel, it was 206,000 MPa.The stud's tensile strength was identified as 496.79 MPa.This numerical demonstration utilized displacement loading, where the axial and shear directions underwent proportional radial loading until reaching a predetermined displacement threshold.The analysis procedure was segmented into 1000 incremental steps, with each step approximating a 0.007 mm increment in the shear direction.Throughout the analysis, the peak incremental step number was capped at 10,000, starting with an initial step size of 0.001.The step size ranged from a minimum of 10 −5 to a maximum of 10 −3 .Various allowable error thresholds ϵ tol , including 1, 10 −2 , 10 −4 , 10 −6 , and 10 −8 , were employed for the calculations.The backward Euler integration method is theoretically capable of delivering precise integration outcomes, with smaller allowable errors yielding more accurate results.Consequently, a scheme adopting an allowable error of ϵ tol = 10 −8 was considered to approach the exact solution most closely.
Figure 13 presents the anticipated load relation curves alongside the load-displacement curves for various integration methods with differing error tolerances ϵ tol .Observations from Figure 13 reveal that within the framework of the backward Euler integration method, setting ϵ tol to one leads to a marked divergence of the calculated results from the precise solution.This discrepancy arises due to the substantial tolerance level failing to adequately amend the elastic forecast in proximity to the yield surface.However, when the allowable error ϵ tol is reduced to 10 −2 or less, the computed outcomes demonstrate high fidelity to the exact solution.To assess the effectiveness of an integration procedure, two statistical indicators can be employed: E F and CPUTIME.The computation of the maximum discrepancy E F in the load for the designated integration process is conducted as detailed below: where F represents the integration load for each respective integration procedure, while F ref serves as the reference load.In this specific case, the computational outcomes derived from the backward Euler method, under an error tolerance set at ϵ tol = 10 −8 , are utilized as the baseline load.CPUTIME indicates the execution time of the integration procedure.When the computational duration for a single macroelement is notably brief, augmenting the number of iterations may contribute to a more precise evaluation of CPUTIME.Developing a graph of CPUTIME and E F against ϵ tol as demonstrated in Figure 14 can provide a visual representation of the integration procedure's accuracy and efficiency.One can observe that CPUTIME remains largely unvaried in response to changes in ϵ tol .Furthermore, the backward Euler integration method, often regarded as theoretically unconditionally stable [48][49][50], renders robust and efficient performance, even when ϵ tol values are minimized.This algorithm's output leans towards convergence without being overly sensitive to predetermined ϵ tol assumptions.However, employing the implicit backward Euler integration within the realm of elastoplastic models necessitates that displacement increments should be maintained in comparatively small sizes to assure convergence.Additionally, executing the implicit algorithm involves greater technical complexity, and introducing kinematic hardening criteria into the elastoplastic model further escalates the numerical computation complexity.

Validation by Experimental Results
To affirm the efficacy of the proposed elastoplastic model for stud connections in accurately predicting load-displacement curves, this study leverages data from push-out tests documented in various research efforts for model validation.Interested readers can refer to the specified publications [41,45,51,52], for detailed experimental data and methodologies.The chosen material properties and geometric parameters mirror those outlined in several referenced studies.An ABAQUS data file was meticulously prepared, and the UEL subroutine, embodying the elastoplastic model, was invoked to simulate the shear force-displacement behavior of the samples.Displayed in Figure 15, the results illustrate that the shear force-displacement curves generated by the elastoplastic model for stud connections closely align with both the experimental and finite element analysis findings, validating the model's capability to accurately predict the load-displacement behavior of monotonically loaded stud connections.Nonetheless, it is important to note that the model falls short of predicting the descending portion of the load-displacement curve.

Conclusions
This research presents a significant advancement in the design and analysis of composite structures by introducing a macroscopic elastoplastic constitutive model for stud connections.The development of this model is rooted in plasticity theory and is informed by empirical data derived from push-out tests as well as finite element simulations.Our results demonstrate that by integrating the combined effects of shear and tensile forces, the model offers a more accurate prediction of the load-displacement behavior of stud connections, which is paramount for ensuring structural integrity and safety.
The comparative analysis with experimental push-out test data revealed that the elastoplastic model offers superior predictive capabilities over traditional models based on nonlinear elastic theory.This highlights the importance of considering elastoplastic effects to more realistically capture the complex behavior of stud connections, especially under conditions that induce plastic deformation.By doing so, the proposed model addresses the limitations of previous models, enabling a more precise analysis of structural stiffness and load-bearing capacity.This precision is crucial for the internal redistribution of forces and the accumulation of plastic damage, which are critical factors in the longevity and safety of composite structures.Moreover, the numerical example implemented in ABAQUS demonstrated the model's stability and robustness, suggesting that the model can be effectively utilized in practical engineering applications.The implementation of the backward Euler method ensures the accuracy of the solution and the consistency of the tangent stiffness necessary for the convergence of the nonlinear finite element analysis.
In conclusion, the macroscopic elastoplastic constitutive model proposed in this study is a substantial contribution to the field of composite structure analysis.It not only enhances the accuracy of predictions concerning the performance of shear connections but also serves as a reliable tool for engineers to design safer and more efficient structures.Future work may extend the application of this model to more complex loading scenarios and explore its integration into broader structural analysis frameworks, thereby reinforcing its value to the engineering community.

Figure 1 .
Figure 1.Positive direction of forces and displacements.

Figure 4 .
Figure 4. Finite element model of push-out test.

Figure 8 .
Figure 8.The definition of the plastic displacement increment angle.

FFFigure 9 .
Figure 9.The interaction between the main program and the subroutine.

Figure 10 .
Figure 10.Possible cases of the loading state.

Figure 14 .
Figure 14.Comparison CPUTIME and E F for different ϵ tol .

Figure 15 .
Verification of the proposed elastoplastic model.(a) First set of data.(b) Second set of data.