Numerical Solution of Bending of the Beam with Given Friction

We are interested in a contact problem for a thin fixed beam with an internal point obstacle with possible rotation and shift depending on a given swivel and sliding friction. This problem belongs to the most basic practical problems in, for instance, the contact mechanics in the sustainable building construction design. The analysis and the practical solution plays a crucial role in the process and cannot be ignored. In this paper, we consider the classical Euler–Bernoulli beam model, which we formulate, analyze, and numerically solve. The objective function of the corresponding optimization problem for finding the coefficients in the finite element basis combines a quadratic function and an additional non-differentiable part with absolute values representing the influence of considered friction. We present two basic algorithms for the solution: the regularized primal solution, where the non-differentiable part is approximated, and the dual formulation. We discuss the disadvantages of the methods on the solution of the academic benchmarks.


Introduction
Mathematical modeling plays a crucial role in the simulations and the development of sustainable civil engineering. This paper aims to present the whole story of the solution process of a given problem, instead of focusing only on one specific step of the solution process. Additionally, we would like to show that in the case of practical computation, in the end, the numerical solver always matters. Indeed, even though the solver has to be chosen in a convenient way referring to the properties of the given problem, the efficiency of the solution crucially depends on the choice of the right numerical approach; actually, there will always be pros and cons, which have to be taken into the consideration.
In the paper, we examine this rule on algorithms for the numerical solution of the beam bending problem of a thin fixed beam with an internal point obstacle with possible rotation and shift depending on a given swivel and sliding friction. The beam is a structural element that primarily resists loads applied laterally to the axis. In the case of a thin beam, we suppose that the diameter of the cross section is much smaller than the length of the beam. The thin fixed beam considered in this paper is the thin beam fixed on both sides.
We suppose that the length of the beam is much bigger than its height. The considered load function applies to small deflections of a beam without considering the effects of shear deformations. Such an example belongs to the most popular benchmarks in engineering practice because of its linearity, and it is well known as the Euler-Bernoulli beam model (see, e.g., in [1], where authors present the derivation of the model). For the nonlinear beam model, which supposes large deformations, we can choose, for example, the Gao beam model, which has been originally proposed in [2]. The model is commonly used both in practice and theory, for example, the extension for contact problems with obstacle [3], friction [4], the extension incorporating nonlocality and surface energy effect [5], or the application to optimal control [6].
From the mathematical point of view, the classical formulation of the problem is represented as a linear differential equation of the fourth-order with corresponding addi-tional conditions. In this case, the weak formulation is defined as an elliptical variational inequality of the second type with non-differentiable functional. For more details see, e.g., the comprehensive books in [7][8][9]. In this paper, we shortly review the formulation in Section 2.1, where we also shortly discuss the existence and uniqueness of the solution in the case of the continuous formulations and the approximated problem.
Convergence analysis and an algebraic formulation is presented in Section 2.4. To solve the problem numerically, we adopt the Finite Element Method (FEM) with Hermite spline functions (see, e.g., the introductory book in [10] or the comprehensive numerical analysis in [11]), see Section 2.6.
Two methods are used to solve the problem in Section 2.7. The first one is the method of regularization, which is based on the approximation of the non-differentiable functional by the differentiable one. The second approach is the dualization of the problem and its reformulation into the so-called Quadratic programming (QP) problem [12]. The corresponding dual problem in our case is a minimization of a strictly convex quadratic function on a feasible set defined by box constraints.
As an example of the presented theory and methods of solution, we consider three benchmarks with increasing difficulty, see Section 2.8. These benchmarks show the basic properties of adopted numerical approaches. The results are discussed in Section 3.
Finally, Section 4 concludes the paper. Please see the section titled "Appendix A" for the list of notations used in the paper.

Problem Formulation
We consider a thin fixed beam of length l > 0 with special internal point obstaclex with possible shift and rotation depending on a given sliding g 1 ≥ 0 and swivel g 2 ≥ 0 friction, see Figure 1. The geometry of the beam is defined by the moment of inertia of the cross section J(x), x ∈ [0, l], and the material of the body is defined by Young's modulus of elasticity E(x). The problem is to find the deflection of the beam u, which is caused by the load function f and limited by the considered type of obstacle.

Classical Formulation of the Problem
To introduce the continuous classical formulation, we assume a beam of length l ∈ R + , load function f ∈ C((0, l)), E, J ∈ C 2 ((0, l)), and internal point obstaclex ∈ (0, l), given sliding friction g 1 ≥ 0 and given swivel friction g 2 ≥ 0. We consider the following problem: find a function u, such that Function u, corresponding to the problem (P TM ), will be called a classical solution of the considered problem.
The differential equation (1) is a mathematical model of the bending of the Euler-Bernoulli beam. The boundary conditions (2) define the clamping of the beam at both ends, relations (3) and (4) represent the conditions of given sliding friction, and relations (5) and (6) are the conditions of the given swivel friction.

Weak and Variational Formulation
As the classical formulation of the problem is too strong and therefore not suitable for direct numerical solution, we formulate the considered problem as a variational one.
Let u be a solution of the problem (P TM ) and v ∈ V. By scalar multiplication of Equation (1) by the function v − u and integration (using Green's theorem-see, e.g., [14]) on the interval (0, l), we obtain From (4) and (6), we have Using (12), (13) and (3), (5), it follows from (11) that Afterwards, we can define a weak formulation of the problem (P TM ) as a problem of finding a function u such that where The sought function u is a weak solution of the problem (P TM ).

Remark 1.
The weak formulation (P TM ) of the problem (P TM ) is formulated as an elliptical variational inequality of the second type with non-differentiable functional q (see, e.g., in [8]).
Obviously, the form a is bilinear and symmetric. To prove its continuity we use Hölder's inequality (see, e.g., in [7]) where c = E L ∞ J L ∞ . From the assumptions (17) and (18) of functions E, J and from v 2 (0) = 0, we have Here, we apply Friedrichs's inequality (see, e.g., [7]) to obtain the V-ellipticity of the form a a(v, v) where c = E 0 J 0 k > 0 and k is the positive constant from Friedrichs's inequality. It is also clear that the functional F is linear. Using Cauchy-Schwarz inequality (see, e.g., in [7]), we can prove its continuity, i.e., F ∈ V * .
Let t ∈ (0, 1). Then, for functional q given by (10), it holds that which means that q is convex. Immediately from the definition of q, we can see that and thus q is proper on V. It is clear that the functional q is continuous. The convexity and continuity of q imply its semi-continuity from below. Therefore, the existence and uniqueness of the solution of problem (P TM ) are guaranteed by the following Theorem 1, the assumptions of which are thus fulfilled.
Theorem 1. Let a : V × V → R be bounded, V-elliptical bilinear form, F ∈ V * and let functional j : V → R be convex, semi-continuous from below and proper on V. Then, there exists an unique solution of For more detail, see, e.g., in [8,15].
As a is bounded; the V-elliptical is bilinear and of symmetric form; functional q is convex, semi-continuous from below, and proper on V; and the problem (P TM ) is equivalent to the variational problem of minimizing quadratic energy functional, i.e., where

Approximation
Let f ∈ L 2 ((0, l)) and E, J ∈ L ∞ ((0, l)) satisfy relations (17) and (18), respectively. Let the form a : V × V → R and the functional F : V → R be given by (15) and (16), respectively. We consider a system of partitions such that points x i , i = 1, . . . , n(h) − 1 are nodal points, i.e., where n(h) + 1 is the number of nodal points of D h and h is the maximum length of intervals T k(h) . Moreover, the point obstaclex ∈ D h for any h. We define the following finite-dimensional subspace V h ⊂ V: Therefore, spaces V h satisfy classical boundary conditions and additionally and lim We approximate the functional q for each h by functional q h : Functional q h is convex, semi-continuous from below and proper in the space V h . We can approximate problem (P TM ) by the sequence of problems of finding u h such that The existence and uniqueness of the solution of (P TM ) h is guaranteed by the Theorem 1.
As the form a is symmetric, the problem (P TM ) h is equivalent to the problem of finding u h such that To approximate the function v ∈ H 2 ((0, l)), we use the Hermitian cubic spline function Additionally, lim Remark 2. For the given function v ∈ H 2 ((0, l)), the Hermitian cubic spline function is given by

Convergence Analysis
To prove the convergence of the approximated solutions u h to the weak solution u of the considered problem, we use the following theorem.
where X ⊂ V such that X = V. Let the system of functionals j h meet the following two conditions: Then, it holds that For more details, see, e.g., in [16].
As the functional q h is convex and semi-continuous from below on V h , it is also weakly semi-continuous from below. In our case, operators r h are operators of the appropriate Hermitian interpolation. From the (28), it follows that Therefore, lim This means that assumptions of Theorem 2 are satisfied, i.e., the sequence u h of solution approximations of the problem (P TM ) h converges to the solution u of the problem (P TM ) in the V-space norm.

The Algebraic Formulation
Let the form a and the functional F be defined by (15) and (16), respectively. We choose To obtain a shape of basis functions, we apply the Hermite's cubic splines. We have  Any function v h ∈ V h can be uniquely written in the form of linear combination of basis functions, i.e., where c i ∈ R are the corresponding coordinates of v h in the basis. It is clear that ϕ i , c i depend on h. If we substitute v h from (33) into the functional J h given by (27), we obtain Then, the problem (P TM ) h can be rewritten in the equivalent form where it is required to find c * such that whereJ T −1 is inverse of T , andĵ is index of the nodal point xˆj =x.
Problem (P TM ) N is the algebraic form of problem (P TM ) h .

Numerical Realization
In this section, we examine two different approaches for the practical numerical solution of the problem. For simplicity, we consider an equidistant grid 0 = x 0 < x 1 < · · · < x n = l, x j = jh, j = 0, . . . , n, where n + 1 is a number of nodal points and h = l/n is a size of intervals. Following the Ritz-Galerkin method, we approximate V by finite-dimensional space Additionally, we assume that there existsĵ such thatx = xˆj, i.e., the point with defined conditions (3)-(6) is one of the used nodal points.

Remark 3.
Note that dim V h = N = 2(n − 1) because the solution in boundary points is given by conditions (2).
The optimization problem (P TM ) N can be written in the form symmetric positive definite (SPD) matrix A ∈ R N,N , and vector b ∈ R N . The solution c * ∈ R N of (38) represents the coordinates (33) of unknown discretized deflection u h in the considered basis.
Although the problem (38) has a simple structure, we are still dealing with a nonlinear optimization problem with a non-differentiable objective function. In our case, we are dealing with absolute values in q N (c) (39). The optimization problems with absolute values in objective functions arise in various applications, and they are already well examined. For example, when one uses L1-regularization (i.e., so-called least absolute shrinkage and selection operator-LASSO [17]) for regularization of linear regression, the problem has a similar form to (38).
Generally, there exist two types of numerical methods to solve the problem-approximation of a non-differentiable term in objective function or dualization of the problem.

Method of Regularization
The idea is to replace (approximate) the non-differentiable term in objective function with a more suitable differentiable one. In our case, we decided to use the piecewise quadratic approximation, i.e., we approximate absolute value ω(z) = |z| bỹ where ε > 0 is a sufficiently small regularization parameter. The example of this approximation is presented in Figure 3. The derivative is given bỹ Using this approximation, the non-differentiable function q N (Φ, c) (39) can be written in formq N,ε (c) = g 1ωε (c 2ĵ−1 ) + g 2ωε (c 2ĵ ).
AS we use piecewise quadratic approximation (40), the derivative is piecewise linear function (41). Additionally, the gradient of quadratic function ψ (39) is linear, and conse-quently the necessary optimality condition for unconstrained optimization problem (38) is given by the system of linear equations The components of gradientq N,ε (c) ∈ R N are equal to zero except for the components corresponding to partial derivatives of c 2ĵ−1 and c 2ĵ . These values depend on conditions from definition of (40) with respect to values c 2ĵ−1 and c 2ĵ . The values are presented in Table 1.
The final system of linear equations can be written in form whereÂ is a matrix A with diagonal elements modified by coefficients of linear terms in ∇q N,ε andb is a vector b with additional constant terms from ∇q N,ε . It is necessary to solve the problem for all possible cases with respect to c 2ĵ−1 , c 2ĵ .
The solution of the problem (P TM ) N (38) is the vector c, for which the corresponding condition on c 2ĵ−1 , c 2ĵ with respect to ε is satisfied.

Dual Problem
The idea of this approach is based on simple observation of the following equivalent form of absolute value: Using this identity, we can rewrite the function q N (Φ, c) (39) into where we denoted Using the separability of variables and the saddle-point property [12], we can write optimization problem (38) in equivalent form .
(47) Function L : R N × 2 → R can be considered as a Lagrange function. From the first Karush-Kuhn-Tucker optimality condition [18], we can derive (note that A is SPD, i.e., non-singular) and substitute into objective function L to get (for more details see in [19]) As the original problem (47) is a maximization problem, we can change the sign of the objective function to obtain Please, notice that the dual problem (49) is a minimization of a strictly convex quadratic function on a feasible set defined by box constraints. Such a problem always has a unique solution [12]. The dimension of the problem is 2, which is much more lower number then the dimension of original primal problem (38). Moreover, the problem belongs to the most basic nonlinear optimization problems and it is solvable by several types of methods, for example, Interior Point, Active set, or Projected Gradient methods [18].

Numerical Experiments
As a demonstration of the presented theory and methods of solution, we consider three numerical benchmarks: with analytical solution (Benchmark 1), with non-trivial load function (Benchmark 2), and with more internal points with given sliding and swivel friction (Benchmark 3).

Benchmark 1
We consider a most simple case: we suppose problem (P TM ) (1)-(6), i.e., optimization problem (P TM ) N (38), with constant E, J, f ∈ R. It can be easily shown that this problem has an analytical solution. In this paper, the solution is used for measuring the absolute error of proposed numerical algorithms. The form of solution is given by with γ = f 24EJ and unknownû 0 = u(x) andû 1 = Du(x). These unknown constants can be computed from the conditions (3), (5), (4), and (6). We derived and simplified the corresponding derivatives and form the system of nonlinear equation and inequalities: We solve the system (51) by solving only the equations (which consist of the elimination of absolute values, dividing by theû 0 andû 1 to obtain systems of linear equations), and afterwards we choose the solution which satisfies all the equations and inequalities (51).
To use the numerical algorithms, we will need to prepare objects in optimization problem (P TM ) N (38). It can be easily shown that in the case of constant E, J, f ∈ R, matrix A is a block tridiagonal Toeplitz matrix and together with vector b, it can be written in a form To provide a specific set of data for benchmark, we examine algorithms on the example of a steel beam (E = 2.15 × 10 11 Nm −2 ) of the length l = 1 m. The point obstacle with given friction is given byx = 8/10 m with friction values g 1 = 10 2 N and g 2 = 5 × 10 4 N. We suppose the load function f = 50,000 N. The beam has a rectangular cross section of height v = 0.02 m and width s = 0.02 m. The analytical solution of this problem is presented in Figure 4. We start our examination with number of used elements as n = 20 and analyze the influence of regularization parameter on absolute error of the solution of regularized problem. We introduce the error measure by whereū is analytical solution and u * ε is numerical solution computed by solving (43) with regularization parameter ε. This formula measures the relative difference between the solutions in used nodal points. The results are presented in Figure 5. In the case of the dual problem, the absolute error decreases naturally with the increase of the discretization density, see Figure 6. The main bottleneck of this approach is the computation of the inverse.

Benchmark 2
In this benchmark, we consider the same problem parameters as in the case of Benchmark 1 except for a shape of cross section and a load function. We consider a solid circular cross section of radius r = 0.01 m and a load function with given parameters f 0 , f 1 ∈ R. Obviously, f is a linear function with f (0) = f 0 and f (l) = f (1) = f 1 .
Computing the corresponding integrals, it is easy to derive the components of a vector of right-hand sides We set f 0 = 10 5 N and solve the problem for various choice of f 1 , see Figure 7 (left). To compare the regularization method and the dual problem method, we solve the problem with f 1 = 1.5 × 10 5 N for increasing the problem dimension n. The computational time is demonstrated in Figure 7 (right). To decrease the small oscillations of the results (especially in the case of short time), we compute all problems 100 times and present the average computational time.

Benchmark 3
In our last benchmark, we increase the number of internal point obstacles with given friction coefficients to 4 equidistantly distributed on [0, 1]. We consider the coefficients of given sliding and swivel friction the same for all points g 1 = 10 4 N and g 2 = 5 × 10 4 N. We consider material constants defined in Benchmark 1 and the solid circular cross section and load function defined in Benchmark 2. The solution provided by the dual problem approach is presented in Figure 8 (left). In the case of 4 obstacle points, the number of corresponding Lagrange multipliers λ is 8, g = [g 1 , g 2 , g 1 , g 2 , g 1 , g 2 , g 1 , g 2 ] T ∈ R 8 and matrix B (45) has 8 rows. The dimension of dual problem (49) is 8; however, the most time-consuming operation remains the inverse of the Hessian matrix, see Figure 8 (right).

Discussion
For the demonstration, we plot the largest condition number from nine systems, which has to be solved, see Figure 5. The exact source of the problem is the term 1/ε added to diagonal elements of the Hessian matrix of corresponding regularized problem, see Table 1 for exact equations. Figure 5 demonstrates the typical problem with regularization: the decrease of regularization parameter decreases the error of solution; however, it consequently increases the condition number of the system matrixÂ. The additional term from the gradient of approximated function (Table 1) is shifting some of the diagonal values of the original SPD matrix A. Using the Gershgorin circle theorem [20], we can state that with decreasing ε, the corresponding centers of Gershgorin discs (which includes eigenvalues) are converging to infinity, keeping the radius of disks (sum of absolute values of non-diagonal elements in the corresponding row) constant. Consequently, some of the eigenvalues are converging to infinity. As other eigenvalues remain unchanged, the condition number of the modified matrixÂ is converging to infinity. The system becomes ill-conditioned.
The speed of convergence of iterative numerical algorithms for solving SPD systems (e.g., Conjugate Gradient method [18]) depends crucially on condition number. On the other hand, if we decide to solve the problem using the direct method (e.g., Gauss elimination method), we hit the finite limits of arithmetical operation. To be more specific, at some point in the process, we have to divide by large number on a diagonal. We can conclude that with decreasing ε, it is harder to solve the system precisely.
On the other hand, the dual algorithm is not using any regularization, and therefore it is not necessary to tune this parameter to obtain the optimal ratio between the condition number of the corresponding matrix and the error of approximation. However, this approach has a big price: the assembly of the dual problem. In dual problem (49), one has to compute the inverse of the Hessian matrix, which typically scales with O(n 3 ). In our benchmarks, the Hessian matrix is sparse with a suitable structure to obtain scaling O(n 2 ), see results from Benchmark 1 presented in Figure 6 (right), as well as solution of Benchmark 3 presented in Figure 8 (right). In comparison to the regularization method, the computation of inverse is at some point the bottleneck of the dual problem approach, and the method is slower than the regularization approach, see Figure 7 (right). It is necessary to mention that there exist techniques to avoid the computation of the inverse, e.g., the factorization or the matrix-free implementation [21]. Additionally, the reconstruction of primal solution from dual solution (48) can be performed as a solution of a system of linear equations instead of the multiplication with a dense inverse matrix. This is the reason why (using our naive implementation) the reconstruction step scales with O(n 2 ) in Figure 8.
It is necessary to highlight that we did not solve Benchmark 3 with the regularization method. For multiple obstacle points, the number of corresponding equations (43) would be equal to the number of all combinations of Table 1 constructed for each obstacle point. This number scales exponentially with the number of obstacle nodes. For instance, in the case of 4 points, we would have to solve 9 4 systems of linear equations (43). On the other hand, the dual problem is a much more efficient approach: increasing the number of obstacle points increases the size of the dual problem linearly.

Conclusions
In this paper, we presented both the theoretical and practical aspects of solving the bending problem of the beam: starting from variational formulation, discussing the discretization by Finite Element Method, and ending with numerical methods for solving the corresponding optimization problem. Our numerical experiments revealed the advantages and disadvantages of the two most basic techniques for dealing with the non-differentiable friction term in objective function: the regularization method and the dual problem approach.
The performance and efficiency of the regularization method depend on the a priori chosen parameter. If this parameter is large, the approximation of the original function is not sufficient, and consequently the solution of the approximated problem can be wrong. If the parameter is small, the corresponding problem is ill-conditioned and hard to solve. Additionally, in the case of the problem with more obstacle points, this approach becomes practically inapplicable.
The method of dual problem transforms the problem onto the problem of Lagrange multipliers. The size of the problem is equal to the number of defined friction conditions in the original problem. However, the transformation requires the computation with the inverse of the Hessian matrix of the original problem.
Nevertheless, the explicit computation of the inverse matrix can be avoided using smart implementation techniques. This will be the topic of our future research. • V * -dual space of V; • (u, v)-scalar product in L 2 ((0, l)).