CC BY
1
2023. 19(5). 491-501 Строительная механика инженерных конструкций и сооружений Structural Mechanics of Engineering Constructions and Buildings
HTTP://JOURNALS.RUDN.RU/STRUCTURAL-MECHANICS
АНАЛИТИЧЕСКИЕ И ЧИСЛЕННЫЕ МЕТОДЫ РАСЧЕТА КОНСТРУКЦИЙ
ANALYTICAL AND NUMERICAL METHODS OF ANALYSIS OF STRUCTURES
DOI 10.22363/1815-5235-2023-19-5-491-501
69o04^ ™ RESEARCH ARTICLE / НАУЧНАЯ СТАТЬЯ
EDN: HSRVMZ
Algorithm for calculating the problem of unilateral frictional contact with an increscent external load parameter
Alexander N. Popov Alexander D. Lovtsov
Pacific National University, Khabarovsk, Russian Federation S SanyaPov@mail.ru
Article history
Received: June 12, 2023 Revised: August 27, 2023 Accepted: September 15, 2023
Conflicts of interest
The authors declare that there is no conflict of interest.
Authors' contribution
Undivided co-authorship.
Abstract. The subject of the study is the contact interaction of deformable elements of building structures. Variational formulations are usually used to solve the problem of modeling of unilateral interaction taking into account friction in the contact zone. An alternative to the popular formulations of discretized problems and iterative methods for their solution is proposed. The problem of contact with friction is expanded in the form of a linear complementarity problem (LCP). To solve the linear complementarity problem, the Lemke method with the introduction of an increasing parameter of external loading is used. The proposed approach solves the degenerated matrix in a finite number of steps, while the dimensionality of the problem is limited to the area of contact. To solve the problem, the initial table of the Lemke method is generated using the contact matrix of stiffness and the contact load vector. The unknowns in the problem are mutual displacements and interaction forces of contacting pairs of points of deformable solids. The proposed approach makes it possible to evaluate the change in working schemes as the parameter of external load increases. The features of the proposed formulation of the problem are shown, the criteria for stopping the stepwise process of solving such problems are considered. Model examples for the proposed algorithm are given. The algorithm has shown its efficiency in application, including for complex model problems. Recommendations on the use of the proposed approach are given.
Keywords: building structures, structural nonlinearity, unilateral links, linear complementarity problem, numerical models, finite element method, increasing load
For citation
Popov A.N., Lovtsov A.D. Algorithm for calculating the problem of unilateral frictional contact with an increscent external load parameter. Structural Mechanics of Engineering Constructions and Buildings.2023;19(5):491-501. http://doi.org/10.22363/1815-5235-2023-19-5-491-501
Alexander N. Popov, Lecturer-researcher at the Higher school of Industrial and Civil Engineering, Pacific State University, Khabarovsk, Russian Federation; ORCID: 0000-0001-6762-5476; E-mail: SanyaPov@mail.ru
Alexander D. Lovtsov, Dr. of Engineering, Professor at the Higher school of Industrial and Civil Engineering, Pacific State University, Khabarovsk, Russian Federation; ORCID: 0000-0001-5050-466X; E-mail: lad@pnu.edu.ru
© Popov A.N., Lovtsov A.D., 2023
(DCS) I ' ls work is licensed under a Creative Commons Attribution 4.0 International License https://creativecommons.Org/licenses/by-nc/4.0/legalcode
Алгоритм расчета задачи одностороннего контакта с трением с нарастающим параметром внешней нагрузки
А.Н. Попов А.Д. Ловцов
Тихоокеанский государственный университет, Хабаровск, Российская Федерация И SanyaPov@mail.ru
Аннотация. Предметом исследования является контактное взаимодействие деформируемых элементов строительных конструкций. Для решения задачи моделирования одностороннего взаимодействия с учетом трения в зоне контакта чаще всего используются вариационные постановки. Предлагается альтернатива популярным постановкам дискретизованных задач и итерационным методам их решения. Задача контакта с трением расширяется в виде линейной задачи дополнительности. Для решения линейной задачи дополнительности применяется метод Лемке с введением нарастающего параметра внешнего нагружения. В предлагаемом подходе решается вырожденная матрица за конечное число шагов, при этом размерность задачи ограничена областью контакта. Для решения задачи формируется начальная таблица метода Лемке с использованием контактной матрицы жесткости и контактного грузового вектора. В качестве неизвестных в задаче выступают взаимные перемещения и усилия взаимодействия контактирующих пар точек, деформируемых тел. Предлагаемый подход позволяет оценить смену рабочих схем по мере роста параметра внешнего воздействия. Показаны особенности предлагаемой постановки задачи, рассмотрены критерии остановки шагового процесса решения таковых задач. Приведены модельные примеры для предлагаемого алгоритма. Алгоритм показал свою эффективность в применении, в том числе и на сложных модельных задачах. Даны рекомендации по использованию предлагаемого подхода.
Ключевые слова: строительные конструкции, конструктивная нелинейность, односторонние связи, линейная задача дополнительности, численные модели, метод конечных элементов, нарастающая нагрузка
Для цитирования
Popov A.N., Lovtsov A.D. Algorithm for calculating the problem of unilateral frictional contact with an increscent external load parameter // Строительная механика инженерных конструкций и сооружений. 2023. Т. 19. № 5. С. 491-501. http://doi.org/ 10.22363/1815-5235-2023-19-5-491-501
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ácek 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
Попов Александр Николаевич, преподаватель-исследователь высшей школы промышленного и гражданского строительства, Тихоокеанский государственный университет, Хабаровск, Российская Федерация; ORCID: 0000-0001-6762-5476; E-mail: SanyaPov@mail.ru
Ловцов Александр Дмитриевич, доктор технических наук, профессор высшей школы промышленного и гражданского строительства, Тихоокеанский государственный университет, Хабаровск, Российская Федерация; ORCID: 0000-0001-5050-466X; E-mail: lad@pnu.edu.ru
История статьи
Поступила в редакцию: 12 июня 2023 г.
Доработана: 27 августа 2023 г.
Принята к публикации: 15 сентября 2023 г.
Заявление о конфликте интересов
Авторы заявляют об отсутствии конфликта интересов.
Вклад авторов
Нераздельное соавторство.
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 [2527], 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 been adopted:
> for forces and displacements normal to the contact surface: compressive force of interaction of points of the contact pair xm- > 0; mutual displacement of points of the contact pair Zni > 0 (Figure 1, a);
> for forces and displacements tangential to the contact surface: if the points of the contact pair are conditionally separated normally to the contact zone, then the interaction forces xTi >0 will create a pair of
forces with a clockwise moment; mutual displacement z>0, if it coincides in direction with x^ >0 (Figure 1, b, c).
Papers [47] and [48] proposed a LCP formulation for frictional contact considering initial gaps:
ZTi > о
xri < 0
zTf > 0
b
Figure 1. Unilateral constraints. The signs' rule for interaction forces x and mutual displacements z
Zni > 0, xzi > 0
zni > 0 xni > 0
a
c
xn
+
xT —
x- _
Rn
Rn
- Rrn: zn
- R - f • R z-
R - f • R "■ТТ J "ЛТ _ z - _
nn ПТ
R ТП + f • Rnn R ТТ + f • RnT
R ТП + f • Rnn - R ТТ + f • Rrn:
RFn - Rnn • n ^Т + f • (RFn - Rnn • - Rrn: • "Л - RFт + f • (RFn - Rnn • П) + Rm • П_
Zn > 0; Xn > 0; zn • Xn = 0;
x+ > 0; x- > 0; z+ > 0; z- > 0; z+T • x+ = 0; z-T • x- = 0,
(1)
+
where: xn, zn are vectors [m*1] of interaction forces and mutual displacements of contact pairs along the normal
to the contact zone; xT = (x+ — xT )/2 is the vector [m*1] of interaction forces of contact pairs along the
tangent to the contact zone; zT = (z + — zT)/2 is the vector [mxl] 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 TT 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; Rpn is a contact load vector (CLV) [m*1] for the links which are normal to the contact zone; RpT is a CLV for the links which are tangential to the contact zone; n 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 + R
F;
z > 0; x > 0; zT • x = 0.
(2)
It was assumed that the external influences are: force influence Rf, kinematic influence Rand temperature influence R. The external influence was divided as follows:
Rf + Ra + R = Rc + p • R,
where: pis the vector increment parameter; Rv,Re are vectors of contact loads from any combination of external influences.
Therefore, assuming that Rv = Rf and R = R + R = 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:
[E R -Rv ]•
= 0,
where E is a diagonal identity matrix [3m x3m].
To initialize the solution process (selection of the leading row), an artificial compression (tightening weight) pc of all unilateral links is introduced, and the initial table takes the form:
[E R -Rv ]• z = pc • e,
_ P_
where: e = [1,1,...,1] is a vector with dimension [3m x1].
(3)
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 basis6 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
EA - 1 N У EI - 1 N-m2
F--1
= 1 N
1 m
0.5 m
0.5 m
Figure 2. Scheme of the beam with unilateral supports
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 steps. To initialize the stepped process, an artificial compression 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 parameterp 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 xt1+ xt2+ xt3+ xt1- xt2- xt3- zn1 zn2 zn3 zt1+ zt2+ zt3+ zt1- zt2- zt3- pRv Rc Min ratio
0 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
1 0 1 0 0 0 0 0 0 0 0.375 -0.75 0.375 0 0 0 0 0 0 -0.688 1
2 0 0 1 0 0 0 0 0 0 -0.1875 0.375 -0.1875 0 0 0 0 0 0 -0.4063 1
3 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 17.7777
4 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
5 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
6 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 17.7777
7 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
8 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 xt1+ xt2+ xt3+ xt1- xt2- xt3- zn1 zn2 zn3 zt1+ zt2+ zt3+ zt1- zt2- zt3- pRv Rc Min ratio
19 10.667 0 0 0 0 0 0 0 0 -2 4 -2 0 0 0 0 0 0 1 10.667
1 7.3333 1 0 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 8.3333
2 4.3333 0 1 0 0 0 0 0 0 -1 2 -1 0 0 0 0 0 0 0 5.3333
3 -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.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
5 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
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
7 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
8 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 xt1+ xt2+ xt3+ xt1- xt2- xt3- zn1 zn2 zn3 zt1+ zt2+ zt3+ zt1- zt2- zt3- pRv Rc Min ratio
9 -5.33333 0 0 0 0 0 0 0 0 1 -2 1 0 0 0 0 0 0 -0.5 -5.33333
1 2 1 0 0 0 0 0 0 0 0 4E-15 -7E-16 0 0 0 0 0 0 -0.5 3
2 -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
3 -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
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
5 -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
6 -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
7 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
8 -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 = 3975 N / 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
Fc=1
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.
Fc=1 F ¿=1
J. pF-10.(6) I
1 ^
Xn2=8.(3) Xn3-5.(3)
z„2=0 zn3-0
Xn3-0 b Zn3-5.(3)
Xn3-5.(3) Zn3=0
Xn3-5.(3) Zn3-0
Xn3-0.5 Zn3-0
Figure 3. Mutual displacements (z) and forces of interaction (x) in a beam
q
„_L - 260 m_J
Fc-
1 1
I
F v
L / n L / n L / n L / n
L - 260 m „
0.13644282
0.045258079
л \ --^-А-тТ
[m] Г —D II II II II
L. stitchina , slippage t
b Q - numerical analytical
Figure 4. The console scheme and results
Figure 5. Scheme of sheet piling in soil
a
F
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 1 = 0.27, soil thickness is t = 1 m .
In this example, the external load parameter increases until a ray solution is obtained at p • Fv = 213518V. 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.
References
1. Migórski S. Optimal Control of History-Dependent Evolution Inclusions with Applications to Frictional Contact. Journal of Optimization Theory and Applications. 2020;185:574-596. https://doi.org/10.1007/s10957-020-01659-0
2. Kikuchi N., Oden J.T. Contact Problems in Elasticity. Philadelphia: Society for Industrial and Applied Mathematics. 1988. https://doi.org/10.1137/L9781611970845
3. Duvaut G., Lions J.-L. Inequalities in Mechanics and Physics. Berlin, Heidelberg: Springer Berlin Heidelberg; 1976. https://doi.org/10.1007/978-3-642-66165-5
4. Glowinski R., Lions J.-L., Trémolieres R. Numerical Analysis of Variational Inequalities. Elsevier Science; 1981. https://doi.org/10.1016/S0168-2024(08)70201-7
5. Kravchuk A.S. The variational method in contact problems. The present state of the problem and trends in its development. Journal of Applied Mathematics and Mechanics. 2009;73(3):351-357. https://doi.org/10.1016/j.jappmathmech. 2009.07.004
6. Bathe K.J., Chaudhary A.A. A solution method for planar and axisymmetric contact problems. International Journal for Numerical Methods in Engineering. 1985;21:65-88. https://doi.org/10.1002/nme.1620210107
7. Klarbring A. A mathematical programming approach to three-dimensional contact problems with friction. Computer Methods in Applied Mechanics and Engineering. 1986;58:175-200. https://doi.org/10.1016/0045-7825(86)90095-2
8. Klarbring A. On discrete and discretized non-linear elastic structures in unilateral contact (stability, uniqueness and variational principles). International Journal of Solids and Structures. 1988;24(5):459-479. https://doi.org/10.1016/0020-7683(88)90002-9
9. Hlavácek J., Haslinger J., Necas J., Lovishek J. Solution of Variational Inequalities in Mechanics. New York: Springer; 1986. https://doi.org/10.1007/978-1-4612-1048-1
10. Cottle R.W., Giannessi F., Lions J.- L. Variational Inequalities and Complementary Problems in Mathematical Physics and Economics. New York: John Wiley; 1979.
11. Eck C., Jarusek J., Krbec M. Unilateral Contact Problems. New York: CRC Press; 2005. https://doi.org/10.1201/ 9781420027365
12. Reclitis G., Raivindran A., Regsdel K. Optimization in the technique. In 2 books, Part. 2. Moscow: Mir Publ.; 1986. Available from: http://padabum.com/d.php?id=39920 (Accessed 12.03.2023)
13. Panagiotopoulos P.D. Inequality Problems in Mechanics and Applications. Convex and Nonconvex Energy Functions. Birkhauser: Boston-Base-Stuttgart; 1985. https://doi.org/10.1007/978-1-4612-5152-1
14. Kravchuk A.S., Neittaanmaki P.J. Variational and Quasi-Variational Inequalities in Mechanics. Dordrecht: Springer Netherlands. 2007. Vol. 147. https://doi.org/10.1007/978-1-4020-6377-0
15. Mitra G., Cottle R.W., Giannessi F., Lions J.L. Variational Inequalities and Complementarity Problems — Theory and Application. The Journal of the Operational Research Society. 1981;32(9):848. https://doi.org/10.2307/2581414
16. Acary V., Brémond M., Huber O. On solving contact problems with coulomb friction: Formulations and numerical comparisons. Advanced Topics in Nonsmooth Dynamics: Transactions of the European Network for Nonsmooth Dynamics. 2018:375-457. https://doi.org/10.1007/978-3-319-75972-2_10
17. Noor M.A., Noor K.I., Rassias M.T. New Trends in General Variational Inequalities. Acta Appl Math. 2020; 170(1):981-1064. https://doi.org/10.1007/S10440-020-00366-2/TABLESZ3
18. Vorovich I.I., Aleksandrov V.M. Contact Mechanics Interaction. Moscow: Fizmatlit Publ.; 2001. (In Russ.) Available From: https://www.rfbr.ru/rffi/ru/books/o_460?FILTER_ID=23@2. (accessed 12.03.2023).
19. Raous M. Art of Modeling in Contact Mechanics. In: The Art of Modeling Mechanical Systems. Springer. 2016;570:203-276. https://doi.org/10.1007/978-3-319-40256-7_4
20. Sofonea M., Matei A. Mathematical Models in Contact Mechanics. Cambridge: Cambridge University Press; 2012. https://doi.org/10.1017/CB09781139104166
21. Souleiman Y., Barboteu M. Numerical Analysis of a Sliding Frictional Contact Problem with Normal Compliance and Unilateral Contact. Open Journal of Modelling and Simulation. 2021;9:391-406. https://doi.org/10.4236/ojmsi. 2021.94025
22. Sofonea M., bin Xiao Y. Well-Posedness of Minimization Problems in Contact. Journal of Optimization Theory and Applications. 2021;188(3):650-672. https://doi.org/10.1007/s10957-020-01801-y
23. Sofonea M., bin Xiao Y. Weak formulations of quasistatic frictional contact problems. Commun Nonlinear Sci Numer Simul. 2021;101:105888. https://doi.org/10.1016/J.CNSNS.2021.105888
24. Sofonea M., Shillor M. Tykhonov Well-Posedness and Convergence Results for Contact Problems with Unilateral Constraints. Technologies. 2021;9(1):1. https://doi.org/10.3390/TECHN0L0GIES9010001
25. Leturcq B., Tallec P. le, Pascal S., Fandeur O., Pacull J. NTFA inspired model order reduction including contact and friction. HAL open science. 2021. Available from: https://hal.archives-ouvertes.fr/hal-03202024 (accessed: 12.03.2023).
26. Antolin P., Buffa A., Fabre M., Fabre M.A. A priori error for unilateral contact problems with Lagrange multipliers and IsoGeometric Analysis. HAL open science. 2018. Available from: https://hal.archives-ouvertes.fr/hal-01710816 (accessed: 12.03.2023).
27. Parisch H. A consistent tangent stiffness matrix for three-dimensional non-linear contact analysis. Int Journal Numerical Methods in Engineering. 1989;28(8):1803-1812. https://doi.org/10.1002/nme.1620280807
28. Sofonea M., Souleiman Y. Analysis of a sliding frictional contact problem with unilateral constraint. Mathematics and Mechanics of Solids. 2017;22(3):324-342. https://doi.org/10.1177/1081286515591304
29. Lu L., Li L., Sofonea M. A generalized penalty method for differential variational-hemivariational inequalities. Acta Mathematica Scientia. 2022;42(1):247-264. https://doi.org/10.1007/S10473-022-0114-Z
30. Arghir M., Benchekroun O. A simplified structural model of bump-type foil bearings based on contact mechanics including gaps and friction. TribolInt. 2019;134:129-144. https://doi.org/10.1016/j.triboint.2019.01.038
31. Peng L., Feng Z.Q., Joli P., Renaud C., Xu W.Y. Bi-potential and co-rotational formulations applied for real time simulation involving friction and large deformation. Computational Mechanics. 2019;64(3):611-623. https://doi.org/ 10.1007/S00466-019-01672-9
32. Lukashevich A.A. Computational modelling of stiffness and strength properties of the contact seam. Magazine of Civil Engineering. 2018;81(5):149-159. https://doi.org/10.18720/MCE.81.15
33. Lukashevich A.A., Rozin L.A. On the decision of contact problems of structural mechanics with unilateral constraints and friction by step-by-step analysis. Magazine of Civil Engineering. 2013;1(36):75-81. https://doi.org/10. 5862/MCE.36.9
34. Gholami F., Nasri M., Kovecses J., Teichmann M. A linear complementarity formulation for contact problems with regularized friction. Mechanism and Machine Theory. 2016;105:568-582. https://doi.org/10.1016/IMECHMACHTHEORY. 2016.07.016
35. Tasora A., Anitescu M. A fast NCP solver for large rigid-body problems with contacts, friction, and joints. Computational Methods in Applied Sciences. 2009;12:45-55. https://doi.org/10.1007/978-1-4020-8829-2_3
36. Zheng H., Li X. Mixed linear complementarity formulation of discontinuous deformation analysis. International Journal of Rock Mechanics and Mining Sciences. 2015;75:23-32. https://doi.org/10.1016/j.ijrmms.2015.01.010
37. Ignatyev A.V., Ignatyev V.A., Onischenko E.V. Analysis of Systems with Unilateral Constraints through the Finite Element Method in the Form of a Classical Mixed Method. Procedia Engineering. 2016;150:1754-1759. https://doi.org/ 10.1016/j.proeng.2016.07.166
38. Morozov N.F., Tovstik P.Y. Bending of a two-layer beam with non-rigid contact between the layers. Journal of Applied Mathematics and Mechanics. 2011 ;75(1):77-84. https://doi.org/10.1016/j .jappmathmech.2011.04.012
39. Berezhnoi D.V., Sagdatullin M.K. Calculation of interaction of deformable designs taking into account friction in the contact zone by finite element method. Contemporary Engineering Sciences. 2015;8(23):1091-1098. https://doi.org/ 10.12988/ces.2015.58237
40. Kudawoo A.D., Abbas M., De-Soza T., Lebon F., Rosu I. Computational Contact Problems: Investigations on the Robustness of Generalized Newton Method, Fixed-Point method and Partial Newton Method. Int J for Computational Methods in Engineering Science and Mechanics. 2018;19(4):268-282. https://doi.org/10.1080/15502287.2018.1502217
41. Averin A.N., Puzakow A.Yu. Design of systems with uniateral constraints. Structural mechanics and structures. 2015;1(10):15-32. (In Russ.) Available from: http://nauteh-journal.ru/files/4403ce61-ea9d-4dcf-abbd-ac90a11aec39 (accessed: 12.04.2023).
42. Fabre M., Pozzolini C., Renard Y. Nitsche-based models for the unilateral contact of plates. ESAIM: Mathematical Modelling and Numerical Analysis. 2021;55:941-967. https://doi.org/10.1051/M2AN/2020063
43. Brogliato B., Kovecses J., Acary V. The contact problem in Lagrangian systems with redundant frictional bilateral and unilateral constraints and singular mass matrix. The all-sticking contacts problem. Multibody System Dynamics. 2020;48(2):151-192. https://doi.org/10.1007/S11044-019-09712-1
44. Streliaiev Yu.M. A nonlinear boundary integral equations method for the solving of quasistatic elastic contact problem with Coulomb friction. Journal of Samara State Technical University, Ser. Physical and Mathematical Sciences. 2016;20(2):306-327. (In Russ.) https://doi.org/10.14498/vsgtu1471
45. Zhiltsov A.V. Modified duality scheme for numerical simulation of the contact between elastic bodies. Mathematical notes of NEFU. 2016;23(4):99-114. (In Russ.) Available from: http://www.mzsvfu.ru/index.php/mz/article/ view/modified-duality-scheme-for-numerical-simulation-of-the-contact-between-elastic-bodies (accessed: 12.04.2023).
46. Landenberger A., El-Zafrany A. Boundary element analysis of elastic contact problems using gap finite elements. Comput Struct. 1999;71(6):651-661. https://doi.org/10.1016/S0045-7949(98)00303-4
47. Popov A., Lovtsov A. Frictional contact problem in building constructions analysis. Magazine of Civil Engineering. 2020;100(8):1-11. https://doi.org/10.18720/MCE.100.!
48. Popov A N., Lovtsov A.D. Frictional contact in the linear complementarity problem with gaps account. Structural mechanics and analysis of constructions. 2021;297(4):36-43. https://doi.org/10.37538/0039-2383.202L4.36.43
49. Lemke C.E. Bimatrix Equilibrium Points and Mathematical Programming. Manage Sci. 1965;11(7):681-689. https://doi.org/10.1287/mnsc.1L7.681
50. Lemke C.E. Some pivot schemes for the linear complementarity problem. Complementarity and Fixed Point Problems. 1978;7:15-35. https://doi.org/10.1007/bfb0120779
