Algorithm for calculating the problem of unilateral frictional contact with an increscent external load parameter

Full Text

1. Introduction One of the important tasks of the strength calculation of building structures is the task of determining the parameters of the stress-strain state (SSS) while changing the parameters of external loading [1]. The object of this study is the contact interaction of deformable building structures under increscent external load. Constructively nonlinear problems have been popular since 1970s in works of Kravchuk, Bathe, Kikuchi, Glowinski [2-6]. Klarbring, Hlaváček and Cottle considered variational formulations [7-10]. A step-by-step algorithm is used in the most popular software systems in case of force incrimination problems for non-linear calculation of building structures. According to this algorithm the loading process is divided by the user into several stages (the method of successive loadings). Iterative methods are used at each stage of loading to determine the increments of the structure’s SSS parameters. It is necessary to solve the problem of contact interaction in these problems at each stage. Such tasks have been popular since the 1980s [11-15] and have maintained their popularity until the present [16-20]. These problems are stated in the form of variational inequalities [21-24], and the following numerical methods were used to solve them: Lagrange multipliers [25- 27], penalty functions [28; 29] and their combinations [30-32]. And other methods using contact finite elements [32-33]; quadratic programming approach [34-36]; finite element methods (Spigot-algorithms) [37-39]; and other [40-46]. The user should specify the following parameters at the stage of a problem modeling: 1. The final value of the external load. 2. The number of loading stages (the value of load increment). 3. Method for solving the contact problem and its parameters. The purpose of this work is to create an algorithm that allows tracking the change of working schemes at the parametric increment in external load. The problem expansion parameter, whose physical meaning was the “tightening weight” in contact pairs, was used in [47] in the algorithm for solving the linear complementarity problem (LPC). In this paper, it was proposed to take the external load growth parameter as the parameter of the problem expansion. This approach enables automating the process of load splitting into stages, within each of them a linear problem can be solved. The following tasks arising from this: 1. Program implementation of the algorithm for solving similar problems. 2. Description of the solution peculiarities. 3. Testing the algorithm. 2. Methods The formulation of the calculation of frictional contact problems proposed below considers a node-to-node contact (contact pair). Let m denote the number of contact pairs. It is assumed that the points in each contact pair are connected by unilateral constraints. The constraint that is normal to the contact zone works only on compression and is enabled when these points are in contact and disabled otherwise. Tangential connection to the contact zone is enabled if the interaction forces are less than the ultimate friction forces and disabled if the interaction forces are equal to the ultimate friction forces. This means that slippage of the contact pair points is not possible when the connection is on, whereas it is possible when the connection is off. The following rule for the use of signs has action of points of the contact pair x >0ni ; mutual displacement of points of the contact pair z >0ni (Figure 1, a); Ø for forces and displacements tangential to the contact surface: if the points of the contact pair zni > 0,xni > 0 б) zτi > 0,xτi > 0 в) zτi >0 a b c Figure 1. Unilateral constraints. The signs’ rule for interaction forces x and mutual displacements z been adopted: zτi > 0 Ø for forces and displacements normal to the contact surface: compressive force of inter are conditionally separated normally to the contact zone, then the interaction forces x >0τi will create a pair of forces with a clockwise moment; mutual displacement z >0τi , if it coincides in direction with x >0τi (Figure 1, b, c). Papers [47] and [48] proposed a LCP formulation for frictional contact considering initial gaps: xn Rnn Rnτ - Rnτ zn x+τ = Rτn + f ⋅Rnn Rττ + f ⋅Rnτ - Rττ - f ⋅Rnτ ⋅ z+τ + x-τ - Rτn + f ⋅Rnn - Rττ + f ⋅Rnτ Rττ - f ⋅Rnτ z-τ RFn - Rnn ⋅ η + RFτ + f ⋅(RFn - Rnn ⋅ η) - Rnτ ⋅ η ; (1) - RFτ + f ⋅(RFn - Rnn ⋅ η) + Rnτ ⋅ η zn x z x 0; x+τ ≥ 0; x-τ ≥ 0; z+τ ≥ 0; z-τ ≥ 0; z+τT ⋅x+τ = 0; z-τT ⋅x-τ = 0, Popov A.N., Lovtsov A.D. Structural Mechanics of Engineering Constructions and Buildings. 2023;19(5):491-501 where: xn,zn are vectors [m×1] of interaction forces and mutual displacements of contact pairs along the normal to the contact zone; xτ = (x+τ -x-τ)/2 is the vector [m×1] of interaction forces of contact pairs along the tangent to the contact zone; zτ = (z+τ - z-τ)/2 is the vector [m×1] of mutual displacements of contact pairs tangentially to the contact zone; Rnn is the contact stiffness matrix (CSM) [m×m] for the constraints in the contact pairs along the normal to the assumed contact zone from a single dislocation of nodes of contact pairs along the normal to the specified contact zone; Rττ is the CSM [m×m] for the links introduced in the contact pairs tangential to the contact zone from the unit dislocation of the contact pairs nodes tangentially to the contact zone; f is a coefficient of friction between the nodes of the contact pair; RFn is a contact load vector (CLV) [m×1] for the links which are normal to the contact zone; RFτ is a CLV for the links which are tangential to the contact zone; η is a vector of mutual initial gaps in the contact zone. Thus, three non-negative variables are required to determine the parameters of the stress-strain state (SSS) in the contact pair: one is responsible for the interaction along the normal and two are responsible for the interaction along the tangent. The system of equations and inequalities (1) can be written in the following reduced form: x = R ⋅z + RF; (2) z ≥ 0; x ≥ 0; zT ⋅x = 0. It was assumed that the external influences are: force influence Rf , kinematic influence RΔand temperature influence Rt . The external influence was divided as follows: Rf +RΔ+Rt =Rc + p⋅Rv, where: pis the vector increment parameter; Rv,Rc are vectors of contact loads from any combination of external influences. Therefore, assuming that Rv =Rf and Rc =RΔ +Rt =0, the parametric incrimination of the force load is modeled. It is proposed to extend the formulation (2) by introducing the parameter p of force load increment: x = R⋅z + p⋅Rv; z ≥ 0; x ≥ 0; zT ⋅x = 0; p ≥ 0. To solve the problem, we use the Lemke method [49], [50]. Then, the initial table of the Lemke method takes the form: x [E R -Rv]⋅ z = 0 , p where is a diagonal identity matrix [3m×3m]. To initialize the solution process (selection of the leading row), an artificial compression (tightening weight) of all unilateral links is introduced, and the initial table takes the form: x [E R -Rv]⋅ z = pc⋅e , (3) p where: e = [1,1,...,1] is a vector with dimension [3m×1]. ANALYTICAL AND NUMERICAL METHODS OF ANALYSIS OF STRUCTURES There are peculiarities of the solution by the Lemke algorithm with an increasing parameter of the external influence. First, it is necessary to select Rv as the leading column. The leading row is chosen by the rule of minimum ratio. Then the standard steps of Lemke algorithm are performed. There are two criteria for stopping the stepped process of the solution: Criterion 1. Suppose that at step k the parameter became more than one: p > 1. The criterion is used if it is necessary to obtain a solution for a given load value (the load at which the contact load vector Rv was formed). In order to get a solution to the linear complementarity problem (LCP) it is necessary to return to the previous step; take p out of the basis and determine the values of basis variables by the formula - p⋅Rv, where p =1. Criterion 2. Suppose that at step k a ray solution is obtained (including the first step of the algorithm) for any non-negative value of pk . In this case, in order to get a solution to LCP it is necessary to take the parameter p out of the basis by choosing the corresponding leading line. Ø If the leading element is not equal to zero, then p should be removed from the basisб and the values of the basis variables should be determined by formula - p⋅Rv, where p> pk . From physical point of view, the result obtained should be interpreted as impossibility to change the working scheme at further load increment. If the ray solution is obtained at pk <1, and the solution is to be obtained for a given load value, then the parameter p=1. Ø If the leading element is equal to zero, then removing the parameter p from the basis leads to the undefined basis variables. In this case, at the current step, the values of the basic variables should be obtained as: - p⋅Rv . From physical point of view, obtained result is interpreted as the transformation of the system into a mechanism in case of a further increment of the load. We obtain the ultimate value of the load parameter corresponding to the transition of the system into a mechanism. 3. Results and discussions Many test problems have been solved to verify the algorithm’s operation. Some of the problems with the description of the algorithm operation are given below. Example 1. The scheme of the problem is a beam on three unilateral supports (Figure 2). The load is a vertical concentrated force in the middle of the right span. The algorithm with an increment external influence parameter is implemented within 2 Figure 2. Scheme of the beam steps. To initialize the stepped process, an artificial compression with unilateral supports is introduced in each assumed contact pair Fc =1. Leading rows and columns are highlighted in gray in all the tables below, the initial table of the problem is shown in Table 1. At the first step (Figure 3, a), the parameter increases to p =10.(6) (Table 2). At this value of the parameter of the external load, the moment of detachment of the left support occurs ( xn1 is eliminated from the basis), which indicates that the interaction force xn1 is equal to zero. At the second step, zn1 should be introduced into the basis. The components of the leading column are negative with the exception of two small positive values. These values are the result of round-off errors, so they are assumed to be equal to zero. Thus, there are no positive components in the leading column, so it means that a ray solution is obtained. In this case, the ray solution can be represented as the impossibility of changing the working scheme with further increase of the load parameter (Figure 3, b). To obtain the solution, the parameter p should be taken out of the basis, and only the variable part of the external influence should be considered: p⋅Fv; Fc =0. The final table is shown in Table 3. For this example, it is possible to obtain the solution of the problem for any value of parameter p. For example: for load p⋅Fv =1, after removing the parameter p from the basis, it is possible to take p=1 and obtain the basis variables respectively -1⋅Rv =1 (Figure 3, c). As can be seen from Table 2, round-off errors can lead to values that are close to zero. In order to stop the algorithm in time, a user-defined parameter of the problem accuracy [47] is introduced. Table 1 Initial table of Lemke algorithm Basis xn1 xn2 xn3 xτ1+ xτ2+ xτ3+ xτ1 xτ2 xτ3 zn1 zn2 zn3 zτ1+ zτ2+ zτ3+ zτ1 zτ2 zτ3 pRv Rc Min ratio 0 1 2 3 4 5 6 7 8 1 0 0 0 0 0 0 0 0 -0.1875 0.375 -0.1875 0 0 0 0 0 0 0.09375 1 10.6666 17.7777 17.7777 0 1 0 0 0 0 0 0 0 0.375 -0.75 0.375 0 0 0 0 0 0 -0.688 1 0 0 1 0 0 0 0 0 0 -0.1875 0.375 -0.1875 0 0 0 0 0 0 -0.4063 1 0 0 0 1 0 0 0 0 0 -0.1125 0.225 -0.1125 -0.5 0.5 0 0.5 -0.5 0 0.05625 1 0 0 0 0 1 0 0 0 0 0.225 -0.45 0.225 0.5 -1 0.5 -0.5 1 -0.5 -0.413 1 0 0 0 0 0 1 0 0 0 -0.1125 0.225 -0.1125 0 0.5 -0.5 0 -0.5 0.5 -0.2438 1 0 0 0 0 0 0 1 0 0 -0.1125 0.225 -0.1125 0.5 -0.5 0 -0.5 0.5 0 0.05625 1 0 0 0 0 0 0 0 1 0 0.225 -0.45 0.225 -0.5 1 -0.5 0.5 -1 0.5 -0.413 1 0 0 0 0 0 0 0 0 1 -0.1125 0.225 -0.1125 0 -0.5 0.5 0 0.5 -0.5 -0.2438 1 Table 2 Table for Step 1 of Lemke algorithm Basis xn1 xn2 xn3 xτ1+ xτ2+ xτ3+ xτ1 xτ2 xτ3 zn1 zn2 zn3 zτ1+ zτ2+ zτ3+ zτ1 zτ2 zτ3 pRv Rc Min ratio 19 1 2 3 4 5 6 7 8 10.667 0 0 0 0 0 0 0 0 -2 4 -2 0 0 0 0 0 0 1 10.667 7.3333 1 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 8.3333 4.3333 0 1 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 5.3333 -0.6 0 0 1 0 0 0 0 0 1E-17 0 0 -0.5 0.5 0 0.5 -0.5 0 0 0.4 4.4 0 0 0 1 0 0 0 0 -0.6 1.2 -0.6 0.5 -1 0.5 -0.5 1 -0.5 0 5.4 2.6 0 0 0 0 1 0 0 0 -1 1.2 -1 0 0.5 -1 0 -1 0.5 0 3.6 -0.6 0 0 0 0 0 1 0 0 1.4E-17 0 0 0.5 -0.5 0 -0.5 0.5 0 0 0.4 4.4 0 0 0 0 0 0 1 0 -0.6 1.2 -0.6 -0.5 1 -0.5 0.5 -1 0.5 0 5.4 2.6 0 0 0 0 0 0 0 1 -1 1.2 -1 0 -1 0.5 0 0.5 -1 0 3.6 Table 3 Table for Step 2 of Lemke algorithm Basis xn1 xn2 xn3 xτ1+ xτ2+ xτ3+ xτ1 xτ2 xτ3 zn1 zn2 zn3 zτ1+ zτ2+ zτ3+ zτ1 zτ2 zτ3 pRv Rc Min ratio 9 1 2 3 4 5 6 7 8 -5.33333 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 -0.5 -5.33333 2 1 0 0 0 0 0 0 0 0 4E-15 -7E-16 0 0 0 0 0 0 -0.5 3 -1 0 1 0 0 0 0 0 0 0 -3.1E-15 1.33E-15 0 0 0 0 0 0 -0.5 -1.8E-15 -0.6 0 0 1 0 0 0 0 0 0 3E-17 -1E-17 -0.5 0.5 0 0.5 -0.5 0 7E-18 0.4 1.2 0 0 0 1 0 0 0 0 0 3E-15 -3E-16 0.5 -1 0.5 -0.5 1 -0.5 -0.3 2.2 -0.6 0 0 0 0 1 0 0 0 0 -2E-15 8E-16 0 0.5 -0.5 0 -0.5 0.5 -0.3 0.4 -0.6 0 0 0 0 0 1 0 0 0 3E-17 -0 0.5 -0.5 0 -0.5 0.5 0 7E-18 0.4 1.2 0 0 0 0 0 0 1 0 0 3E-15 -0 -0.5 1 -0.5 0.5 -1 0.5 -0.3 2.2 -0.6 0 0 0 0 0 0 0 1 0 -2E-15 8E-16 0 -0.5 0.5 0 0.5 -0.5 -0.3 0.4 Example 2. An analytical solution of the following frictional contact problem at known friction limit forces was obtained in [1]. A cantilevered beam of length L=260 m was considered, which was placed on a rigid base (Figure 4, a). The load is a constant pressing vertical uniformly distributed load q =3975N /m, the longitudinal force at the right end Fv =243750N was considered as a load with increasing parameter. The analytical solution is compared with the numerical solution obtained by the proposed method. The foundation is modeled by a set of discrete rigid supports. A plane frame finite element (FE) with three degrees of freedom at a node is used to model the console. The cantilevered beam is divided into n=10 elements. The FE nodes contact the supports according to the Coulomb friction scheme. The concentrated vertical pressing force in the node is Fc = q ⋅ L =114833.33 N . n The friction coefficient is assumed to be f = 0.3. The results for the longitudinal displacements of the beam are shown in Figure 4, b. The dependence of the error on the number of accepted elements is shown in Table 4. The error is calculated by the formula: uerr = abs(analytical - numerical)/analytical . Due to round-off errors, there is no clear correlation. Table 4 Error calculation n 5 10 20 40 80 160 U, error, % 0.016 0.085 0.036 0.91 0.021 0.84 Example 3. The problem of plane deformation of a sheet pile wall in soil with an underlying layer of rocky soil is considered (Figure 5). The sheet pile wall interacts with the soil on its two sides according to the Coulomb friction scheme. The sheet pile wall and the soil are conventionally separated in Figure 5. A horizontal concentrated force at the top of the sheet pile wall is taken as a variable load, which is affected by the increasing external load parameter; the dead weight of the soil is not considered. Horizontal displacements for the soil are forbidden on the sides, vertical and horizontal displacements are forbidden at the base, and the pile has a hinge immobile support at the base. A frame element of plane problem with 3 degrees of freedom at a node is used to model the sheet pile wall, and a 4 node element of plane problem of elasticity theory is used for soil. The coefficient of friction between steel and soil is assumed to be f = 0.4. c) Figure 3. Mutual displacements (z) and forces of interaction (x) in a beam Figure 4. The console scheme and results Figure 5. Scheme of sheet piling in soil 497 ANALYTICAL AND NUMERICAL METHODS OF ANALYSIS OF STRUCTURES Popov A.N., Lovtsov A.D. Structural Mechanics of Engineering Constructions and Buildings. 2023;19(5):491-501 The area of cross section of one meter of sheet pile is Ap =1.831⋅10-2 m2 , modulus of elasticity is Ep = 2 ⋅1011 Pa , moment of inertia Ix = 1.016 ⋅10-5 m4 . For soil modulus of elasticity is E = 4.5⋅107 Pa , Poisson's ratio is μ=0.27, soil thickness is t =1m . In this example, the external load parameter increases until a ray solution is obtained at p⋅Fv =213518N. The zones of sheet pile detachment from soil appear at the top of the sheet to the left and at the bottom of the sheet pile to the right. Zones of contact appear at the top of the sheet pile to the right and at the bottom of the sheet pile to the left. In this case, the adhesion zone occurs only on the left side in two nodes. On the right, the soil slides along the sheet pile (Figure 6). Figure 6. Results for Problem 3. Interaction forces x, mutual displacements z It should be noted that according to the results of solving the testing problems, the following feature of the algorithm was revealed. During the step-by-step process of the Lemke algorithm, the problem of comparison with zero arises. The occurrence of small values is due either to the “physics” of the problem (small load increment leading to a change in the working scheme of the structure) or to round-off errors. This leads to the problem of finding a criterion for the difference of these small values from zero. For this purpose, a single artificial parameter for the accuracy of stopping the step-by-step process was introduced. It determines how much the obtained value of the ultimate desired external load p⋅Fv will differ from the exact value within the framework of the discretized problem. The necessity of comparison with zero appears, as a rule: 1) at the last stage if several variables, including the parameter p, tend to leave the basis simultaneously; 2) in the case if the values are close to zero in the leading column or close to zero and negative in the load column. In the first case, one should act according to Criterion 2 for stopping the step-by-step process. In the second case, small values are interpreted as the result of round-off errors and should be assumed to be zero. It has been experimentally determined that the optimum range for the value of the parameter is from 10-14 to 10-4 in the most difficult cases. This describes an absolute error in external load increment parameter. The examples presented in the paper have been selected, among other reasons, to show the effect of round-off errors on the interpretation of the algorithm’s solution results. Thus, in Example 1 (Table 2) small values appeared in the leading column, and in Example 2 (Table 4) the tendency of the numerical solution to the analytical one is not monotonous. due to rounding errors. The value of the external load that gives the ray solution in Example 3 due to the accumulation of round-off errors, it was necessary to decrease the parameter for stopping the step process to 10-4, which did not affect the accuracy of the solution. 4. Conclusions An algorithm for tracking the change of working schemes at parametric increment of external load for structurally nonlinear contact problems with friction has been developed. The algorithm has shown its effectiveness in solving problems with large contact interaction forces. The physical meaning of the algorithm is a sequential change of working schemes (differing one from another by switching of unilateral constraints) at parametric increment of force load. This enables to automate the process of load dividing into stages, within each of which a linear problem is solved. The use of the proposed approach makes it possible to fulfill strictly the condition of mutual non-penetration of contacting bodies. However, if there is a frequent change of working schemes with a small increase in the parameter of external influence, then this leads to the accumulation of round-off errors and to the complication of determining of the criterion for stopping the step process. The algorithm shows good results for problems with a small contact area and large interaction forces in the assumed contact area. The accuracy of calculating the results remains high enough even in difficult conditions for the algorithm. In the process of the work the following tasks have been fulfilled: 1. A Python program has been written that implements the Lemke algorithm with an increment parameter of external influence. 2. A number of features of the algorithm solution have been described, i.e., the beginning of the step process of the solution, its completion and interpretation. 3. The process of solution has been shown and described for a number of testing problems. The peculiarities of the algorithm operation have been identified and shown in examples.

About the authors

Alexander N. Popov

Pacific National University

Author for correspondence.
Email: SanyaPov@mail.ru
ORCID iD: 0000-0001-6762-5476

Lecturer-researcher at the Higher school of Industrial and Civil Engineering

Khabarovsk, Russian Federation

Alexander D. Lovtsov

Pacific National University

Email: lad@pnu.edu.ru
ORCID iD: 0000-0001-5050-466X

Dr. of Engineering, Professor at the Higher school of Industrial and Civil Engineering

Khabarovsk, Russian Federation

References