- Email: [email protected]
Improved algorithm for the detection of bifurcation points in nonlinear finite element problems
Computers and Structures 191 (2017) 1–11
Contents lists available at ScienceDirect
Computers and Structures journal homepage: www.elsevier.com/locate/compstruc
Improved algorithm for the detection of bifurcation points in nonlinear finite element problems S. Léger ⇑, J. Haché, S. Traoré Département de mathématiques et de statistique, Pavillon Rémi-Rossignol, 18 avenue Antonine-Maillet, Université de Moncton, Moncton E1A 3E9, Canada
a r t i c l e
i n f o
Article history: Received 14 October 2016 Accepted 6 June 2017
Keywords: Finite element method Moore-Penrose continuation method Schur complement Mesh adaptation Algorithm Bifurcations
a b s t r a c t Dealing with bifurcation points when solving large deformation finite element problems is not an easy task. Near such points, the Jacobian matrix becomes singular and the problem becomes difficult to solve numerically. In these situations, increasing heuristically the loading parameter during the simulation in order to follow the solution branch is not an option as this approach usually results in the divergence of the process. Efficient numerical techniques capable of handling the presence of bifurcation points are therefore necessary and continuation methods have proved to be powerful tools when dealing with these kind of issues. In Léger et al. (2015), a new implementation technique based on a Schur complement approach for the Moore-Penrose continuation method, which facilitates the detection of bifurcation points and enables branch following, was presented. This method has proved to perform well in most situations; however, in others (i.e. when mesh adaptation is added to the algorithm), some difficulties appear. In this paper, we therefore present an improved approach, which is much more robust, for the detection of bifurcation points in nonlinear finite element problems. Numerical examples will be presented to show the efficiency of the new approach. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction In recent years, the use of numerical methods to solve industrial problems is constantly increasing. With the complexity of these problems constantly increasing as well, developing robust, efficient and accurate numerical methods is a priority. This is especially true for large deformation hyperelastic problems, where convergence problems are frequently encountered. These problems can be related to different issues, but in most cases they are due to distorted elements that appear in the mesh or to the presence of bifurcation points (where the Jacobian matrix is singular). Bifurcation points, which correspond to eigenvalues of the system, are quantities of interest [26,21,27]. Detecting such points and following the different solution branches passed the critical loads associated with the bifurcation points are thus important and useful to better understand the physical properties of the problem we are solving. The analysis of structural instabilities [6,28], which includes the detection of bifurcation points, is generally based on a numerical continuation method (see [9,4,11,15,23,1,22,24]). In a finite element context, where large systems are frequently encountered, simply looking at the sign of the determinant of the
⇑ Corresponding author. E-mail address: [email protected] (S. Léger). http://dx.doi.org/10.1016/j.compstruc.2017.06.002 0045-7949/Ó 2017 Elsevier Ltd. All rights reserved.
Jacobian matrix for the detection of bifurcation points in not an option as the computation of these determinants end up with machine overflows. A new implementation technique for the Moore-Penrose continuation method, which takes advantage of information already available and facilitates the detection of such points in a finite element context, was therefore introduced in [17]. To validate this new algorithm, the classical elastic beam buckling problem was considered as the bifurcation modes for this problem are well known (see Le Tallec and Vidrascu [16], Ikeda et al. [14]). The finite strain incompressible elasticity problem presented in [2,3] was also used for the validation process. Not only were we able to detect the bifurcation point presented in those papers, but we were also able to detect many more bifurcation points which will be presented in this paper. We note that, at that point, the validation process was done using fairly regular meshes. To solve large deformation problems with the finite element method, two formulations are frequently used in practice: the total Lagrangian formulation [13] which refers to the initial configuration and the updated Lagrangian formulation [5] which refers to the most recently calculated configuration. These two formulations are mathematically equivalent, but the updated Lagrangian formulation tends to perform better in practical applications, as was shown in [18]. It was also shown in [18,19] that remeshing the initial configuration, which does not necessitate any transfer of vari-
2
S. Léger et al. / Computers and Structures 191 (2017) 1–11
ables (meaning that no loss of information results from this), also enhances the performance of the updated Lagrangian method. However, when trying to incorporate the detection of bifurcation points in our complete updated Lagrangian algorithm, we noticed that remeshing had an effect on the detection of bifurcation points. Investigation led us to a better understanding of the phenomena observed and modifications were made to the bifurcation point detection algorithm in order to obtain a much more robust numerical method, which is the object of this paper. The paper is organized as follows. Section 2 describes the implementation technique used for the Moore-Penrose continuation method and explains how this approach facilitates the detection of bifurcation points in a finite element context. Section 3 compares the results obtained when using fairly regular meshes versus adapted meshes. Section 4 describes the improved algorithm for the detection of bifurcation points while Section 5 is devoted to the validation of this improved strategy. 2. Moore-Penrose continuation method When solving large deformation problems, a loading parameter k, corresponding to either external forces or prescribed displacements or both, implicitly drives the deformation. This parameter is typically increased gradually until the desired loading, kmax , is attained. However, increasing heuristically the loading parameter while using a standard Newton method to solve the nonlinear system of equations is highly inefficient as this approach will usually result in the divergence of the algorithm in the neighbourhood of limit points or bifurcation points. A more efficient approach is to use a continuation method where the loading parameter k is explicitly introduced in the system of nonlinear equations, which can be expressed as:
FðxÞ ¼ Fðu; kÞ ¼ 0
ð1Þ
with F : RNþ1 ! RN a smooth function. The vector of unknowns, x, therefore consists of the N degrees of freedom denoted as u plus the loading parameter. The Moore-Penrose continuation method, which was shown to be very robust in the case of large deformation problems (see [18]), is a predictor-corrector method and can be summarized as follows (see [10,17]). Starting from a known point xðiÞ 2 RNþ1 on the solution curve, and given a tangent vector v ðiÞ at that point, the next point xðiþ1Þ 2 RNþ1 on the solution curve as well as the tangent vector at that point v ðiþ1Þ can be obtained using the following algorithm: X 0 ¼ xðiÞ þ hv
AðX
k
Þ dkx
> V k dkx
¼0
AðX k ÞT k ¼ AðX k ÞV k >
Vk Tk ¼ 0 3. Update the solution vector:
X kþ1 ¼ X k dkx 4. Update the tangent vector: k
Z ¼V T
k
6. If kFðX k Þk 6 eF and kX kþ1 X k k 6 ex , convergence attained:
xðiþ1Þ ¼ X kþ1 ;
v ðiþ1Þ ¼ V kþ1
We note that X 0 is the prediction point obtained by making a prediction step of length h in the tangential direction, AðX k Þ ¼ ½F 0u ðX k Þ F 0k ðX k Þ is a rectangular matrix of dimension h i> n ðn þ 1Þ, dkx ¼ dku dkk is a correction vector of dimension ðn þ 1Þ 1; kmax represents the maximum number of iterations allowed while eF and ex are the desired tolerances on F and x respectively. As for T k , it represents a correction on the tangent vector. More details can be found in [17]. As can be seen, to obtain the next point on the solution curve, two systems need to be solved at each iteration. In matrix form, these systems are given by:
F 0u ðX k Þ F 0k ðX k Þ V ku
>
F 0u ðX k Þ F 0k ðX k Þ >
dku
!
¼
dkk
V kk
and
V ku
!
!
T ku
!
T kk
V kk
¼
FðX k Þ 0
!
AðX k ÞV k
ð4Þ
!
0
ð5Þ
where the vectors V k ; T k and dkx were decomposed into their comh i> > ponents in the u and k directions (i.e. V k ¼ V ku V kk ; T k ¼ h i> h i> T ku T kk and dkx ¼ dku dkk ). In a finite element context, constructing and solving these two systems can be very costly. A new algorithm, based on a Schur complement approach, was therefore introduced in [17]. This new approach not only solves the systems by using as much as possible information that is already available, but it can also be incorporated numerically in a finite element code by simply adding a postcondition to the Newton method solver. The following algorithm summarizes the implementation technique of the Moore-Penrose continuation method when using this approach. X 0 ¼ xðiÞ þ hv
ðiÞ
V ¼v For k ¼ 0; 1; 2; . . . ; kmax 0
ðiÞ
(b) Solve for W ku :
ð2Þ
F 0u ðX k Þ W ku ¼ F 0k ðX k Þ (c) Calculate the Schur complement bk :
2. Solve the linear system:
(
Z kZk
F 0u ðX k Þ Dku;k ¼ FðX k Þ
k
¼ FðX Þ
V kþ1 ¼
1. (a) Solve to obtain the standard Newton correction Dku;k (for a fixed value of k):
ðiÞ
V 0 ¼ v ðiÞ For k ¼ 0; 1; 2; . . . ; kmax 1. Solve the linear system:
(
5. Normalization of the tangent vector:
bk ¼ V kk V ku W ku ð3Þ
(d) Calculate the Moore-Penrose correction for the load parameter k:
dkk ¼
1 bk
V ku Dku;k
(e) Apply the Moore-Penrose correction on the standard Newton correction:
3
S. Léger et al. / Computers and Structures 191 (2017) 1–11
dku ¼ Dku;k W ku dkk 2. (a) Calculate the correction for the k component of the tangent vector:
1
T kk ¼
bk
½V ku V ku þ V ku W ku V kk
V ku
¼
þ
W ku
½V kk
3. Numerical results: regular meshes vs adapted meshes Let us consider the finite strain incompressible elasticity problem presented in [2,3], where a square of initial configuration
(b) Calculate the correction for the u component of the tangent vector:
T ku
can be used to follow a specific bifurcation branch (by using it as a prediction direction instead of the usual tangent vector).
T kk
3. Update the solution vector:
X kþ1 ¼ X k dkx 4. Update the tangent vector:
Z ¼ Vk Tk 5. Normalization of the tangent vector:
X0 ¼ ½1; 1 ½1; 1, with coordinates x and y and meshed as in Fig. 1, is considered. For this problem, the top edge of the configuration is traction free, while all other edges are fixed. The body is subjected to a vertical uniform body force F, the incompressible elastic material is described by a Mooney-Rivlin model with parameters k ¼ 40; 000; 000; c1 ¼ 20 and c2 ¼ 0 and P2 P1 elements are used for the computation. By using a standard Newton approach to solve this problem, the algorithm converges to the trivial solution, given by u ¼ 0 and p ¼ Fð1 yÞ, until convergence problems are encountered at a load of approximately F ¼ 265, which corresponds to the critical value detected in [3]. In [17], the algorithm presented in Section 2
Z kZk
V kþ1 ¼
6. If kFðX k Þk 6 eF and kX kþ1 X k k 6 ex , convergence attained:
xðiþ1Þ ¼ X kþ1 ;
v ðiþ1Þ ¼ V kþ1
It was also shown in [17] that this new approach facilitates the detection of bifurcation points. In fact, according to [12], by defining sðxÞ via the system:
JðxÞ B C> 0
YðxÞ sðxÞ
¼
0 ; 1
ð6Þ
where JðxÞ is the Moore–Penrose matrix given by:
JðxÞ ¼
; >
A
v
and where B and C are chosen in such a way that B R ImðJð xÞÞ and C R Im J ð xÞ> (with x the bifurcation point), then sðxÞ will change sign at the bifurcation point x ¼ x. By also solving this additional system using a Schur complement approach, the solution is easily obtained and is given, at each iteration k, by:
Fig. 1. Finite strain incompressible elasticity problem: initial mesh.
bk ¼ V kk V ku W ku (Schur complement) Gkk ¼ b1k ðBk V ku P ku Þ Gku ¼ P ku Gkk W ku h i1 sk ¼ C k Gkk þ C u Gku
0.5
Y kk ¼ sk Gkk
0.3
Gku
0.2 1
0.1
1
where P ku ¼ F 0u ðX k Þ Bu and W ku ¼ F 0u ðX k Þ F 0k ðX k Þ and where the vectors B; C and Y have been decomposed into their components in the u and k directions respectively (i.e. B ¼ ½Bu Bk > , C > ¼ ½C > u C k and Y ¼ ½Y u Y k > ).
As bk and W ku are already known, the detection of bifurcation points can be done by simply solving one more linear system with the matrix F 0u ðX k Þ (to calculate Pku ) and by computing a few scalar products. We note that F 0u ðX k Þ is nothing but the classical Jacobian matrix in a classical Newton method. The detection of bifurcation points therefore only requires adding a few lines of code in the postcondition to the Newton method solver. We also note that while s can be used to monitor bifurcation points, the vector Y
L2
¼ s
k
u
Y ku
0.4
0 -0.1 -0.2 -0.3 -0.4 -0.5 0
200
400
600
800
1000
1200
1400
1600
1800
F Fig. 2. Finite strain incompressible elasticity problem: detected bifurcation points.
4
S. Léger et al. / Computers and Structures 191 (2017) 1–11
Fig. 3. Finite strain incompressible elasticity problem: first six bifurcation modes.
3.5 (a)
3 2.5
(b)
u
L2
2 1.5
(c)
1
(d)
0.5 (e) (f)
0 -0.5 0
200
400
600
800
F
1000
1200
1400
1600
1800
Fig. 4. Finite strain incompressible elasticity problem: bifurcation diagram of the first six bifurcation modes.
was used and a bifurcation point was detected at an approximate value of F ¼ 264, which explains the convergence problems of the standard Newton approach. It was also shown in [17] that another bifurcation point exists for this problem at an approximate value of F ¼ 370. Detecting this bifurcation point was made possi-
Fig. 5. Finite strain incompressible elasticity problem: deformed meshes for the 2nd antisymmetric and symmetric modes of bifurcation.
ble by staying on the fundamental branch passed the first bifurcation point (instead of following the bifurcation branch). At that time, a total body force of F ¼ 500 had been imposed for the computation, leading us to believe that only two bifurcation points existed for this problem. However, by simulating this problem for larger values of F, many more bifurcation points were detected by using our proposed algorithm, as can be seen in Fig. 2. We note that the curve
S. Léger et al. / Computers and Structures 191 (2017) 1–11
Fig. 6. Finite strain incompressible elasticity problem: initial meshes used for the comparison.
that passes through the detected bifurcation points (indicated by the symbol () in the figure) represents the fundamental branch. The bifurcation modes associated with these points, which were validated by looking at the eigenvalues, can also be followed by using the vector Y as described in Section 2. Fig. 3 illustrates the first six bifurcation modes detected for the finite strain incom-
5
pressible elasticity problem and also indicates the approximate load at which these different bifurcations occur. As can be seen, the modes associated with the different bifurcation points alternate between antisymmetric and symmetric. A bifurcation diagram of these six bifurcation modes is shown in Fig. 4, while Fig. 5 illustrates the robustness of our Moore-Penrose continuation method by showing the high level of deformation attained for both the 2nd antisymmetric and symmetric modes of bifurcations. We note that similar results were also obtained for the other modes. Now, as remeshing the initial configuration in an updated Lagrangian context enhances the performance of the numerical method, let us remesh the configuration of Fig. 1 by using the fully optimal mesh adaptation method developed by Bois et al. [7,8] in order to verify that the proposed algorithm still detects the bifurcation points when remeshing is added to the algorithm. We note that for this example, the error indicator not only takes into consideration the displacement. It is also based on the pressure due to the fact that the trivial solution for the displacement is null. By taking into consideration both the displacement and the pressure, we ensure that the error indicator is able to detect the variations in the solution for both eventualities (i.e. when the calculated point lies on the fundamental curve and only the pressure varies as well as when the calculated point lies on a bifurcation branch in
Fig. 7. Finite strain incompressible elasticity problem: calculated values of s near the second bifurcation point at an approximate load of F ¼ 370 (comparison of the values obtained using a regular mesh versus an adapted mesh).
6
S. Léger et al. / Computers and Structures 191 (2017) 1–11
which case the displacement is nonzero). Since the initial mesh is adapted based on the error indicator calculated at the end of the first loading step, we note that both eventualities are feasible depending on the size of the loading step (as well as the step size chosen for the Moore-Penrose continuation method). Adding the remeshing of the initial configuration in our algorithm led us to notice that the detection algorithm did not always work properly. To better understand the reason behind this, let us take a closer look at the values of s calculated during the simulations. For the comparison, we will solve the problem with two different meshes. In one situation, the algorithm will detect the bifurcation point, while in the other, the algorithm will not be able to detect the bifurcation point. Fig. 6 illustrates the two different meshes used for the comparison, while Fig. 7 shows the calculated values of s near the second bifurcation point which is located at an approximate load of F ¼ 370. Before comparing the values of s, let us first remind the reader that while following the solution curve, vector G can become orthogonal to vector C which produces a singularity in the computation of s and indicates that the matrix MðxðsÞÞ, defined by
MðxðsÞÞ ¼
JðxðsÞÞ B C
>
0
ð7Þ
is singular for some value of s. Numerically, the algorithm will see this as a change of sign even though s is not defined there. In [17], we had therefore implemented a second verification test in order to eliminate these false bifurcation points. However, we can easily imagine that difficulties can appear if ever a real bifurcation point is located just after a discontinuity of s (or vice versa). In fact, in that specific case, two changes of sign would happen in a very small interval, meaning that the algorithm could very well not see these two changes of sign, all depending on which points are calculated along the solution curve. By looking at Fig. 7, we can see that the discontinuity of s, which occurred at a load of approximately F ¼ 317 in the case of the regular mesh, now occurs at a load of approximately F ¼ 365 in the case of the adapted mesh, meaning that the discontinuity is now much closer to the real bifurcation point and that the scenario described above could easily happen. To confirm that this is the reason why the bifurcation point has not been detected by the algorithm, we have approximated the values of s at different locations near the bifurcation point. Fig. 8 illustrates these values (indicated by the symbol ()) and confirms our hypothesis. The two consecutive converged points near the bifurcation point are respectively located just before the discontinuity of s and just after the real bifurcation, meaning that both values are positive and that no change of sign is detected. This phenomenon seems to happen more frequently with adapted meshes, or finer meshes, as the discontinuity always seems to take place much closer to the real bifurcation point in those cases. As the algorithm cannot dictate which points are calculated by the continuation method along the solution curve, changes are needed in order to ensure that the bifurcation points are detected no matter which mesh is used initially. This improved and much more robust algorithm is described in the next section. 4. Improved algorithm for the detection of bifurcation points Not many papers deal with the detection of bifurcation points in large deformation problems. Ligursky and Renard [20] have however recently published a paper in which they propose a criterion for detecting bifurcation points as well as a procedure for its realization in a numerical continuation context designed for large problems. Their approach is similar to the one described in Section 2, with the exception that they do not seem to use a Schur
complement approach to solve the different systems. Another difference in their approach is that between each consecutive converged points, they approximate matrix J by using the linear combination:
J :¼ ð1 aÞJðxðiÞ Þ þ aJðxðiþ1Þ Þ where a varies from 0 to 1, xðiÞ is the last converged point and xðiþ1Þ is the current converged point, and then solve the system
J
B
C>
0
Y
s
¼
0 ; 1
ð8Þ
for each a in order to approximate the values of s in between the points. They do not specify that this is to help with the detection of bifurcation points when two changes of sign happen in small intervals, but we can certainly assume that it is, as this procedure gives a better idea of the behaviour of s in between the converged points (as can be seen in Fig. 8 where this procedure was used to approximate the values of s). Their algorithm is described in Algorithm 1. Algorithm 1. Ligursky and Renard’s algorithm [20] for the detection of bifurcation points Require: xðiÞ ; xðiþ1Þ ; v ðiÞ ; v ðiþ1Þ ; dmax > dmin > 0; dinc > 1 > ddec > 0 Step 1: Set nch :¼ 0; a :¼ 0; s0 :¼ 106 ; choose B and C randomly. Step 2: Set
s1 :¼ 106 ; d :¼ dmin , and
J :¼ ð1 aÞJðxðiÞ Þ þ aJðxðiþ1Þ Þ;
M :¼
J C>
B 0
Step 3: Solve
M
Y
s
¼
0 1
Step 4: If ss0 < 0 and js0 j < js1 j, set nch :¼ nch þ 1. Step 5: If js s0 j is large, set d :¼ maxfddec d; dmin g, otherwise if js s0 j is small, set d :¼ minfdinc d; dmax g. Step 6: If a < 1, set a :¼ minfa þ d; 1g; s1 :¼ s0 ; s0 :¼ s and return to Step 2. Step 7: If nch is odd, a bifurcation has been detected. By incorporating this algorithm in our procedure, the bifurcation points are more easily detected. In fact, by using this modified approach (i.e. approximating the values of s in between the converged points), the algorithm now detects the bifurcation point located at F ¼ 370 in the case of the adapted mesh. However, this modification does not completely solve the problem as can be seen in Fig. 9(a) where a bifurcation point is not detected even with the use of the approximated values. In this specific case, a discontinuity of s seems to occur in the interval of interest, but the location of the intermediate points does not enable us to detect the change of sign. By taking smaller values of dmax and dmin , we are able to better approximate the behaviour of s and detect the bifurcation point as can be seen in Fig. 9(b). This however has a fairly high cost associated with it. Instead of taking 0.55 s (as in the first set of parameters) in this specific interval, it takes 3.58 s. As this has to be done between each consecutive converged points, reducing the values of the parameters for the entire simulation would be very costly and
S. Léger et al. / Computers and Structures 191 (2017) 1–11
7
Fig. 8. Finite strain incompressible elasticity problem: approximation of the values of s near the real bifurcation point in the case of the adapted mesh. The approximated values are indicated by the symbol (), while the converged points are indicated by the symbol ().
is not recommended. Furthermore, the number of intermediate points necessary to detect the bifurcation point depends on a number of different aspects (i.e. mesh, problem we are solving, distance between the two converged points, etc.). It is therefore not possible to find optimal values that would work in all cases, unless we take extremely small values for the parameters, which again would result in higher computational times. For example, by using the same values of parameters as in [20], which are extremely small, it takes 514.8 s to detect the bifurcation point in this specific interval. As bifurcations usually only occur a few times during a simulation, there is no need to calculate many intermediate points between each consecutive converged points. However, we do need to calculate a few intermediate points each time as the bifurcation points can easily be hidden when discontinuities of s occur (as was shown in Fig. 8). In all other cases, the intermediate points usually follow a smooth trajectory and lead to the same result as the original algorithm for the detection of bifurcation points, as can be seen in Fig. 10, where in this case, the values of s remain strictly increasing even when more points are calculated. Taking smaller values for the parameters are unnecessary in those cases. Our proposed strategy therefore consists of using larger values of parameters throughout the simulation as to get an idea of the behaviour between the points without generating a high computational cost and then only refining these parameters when needed. To determine if further investigation is needed, we study the shape of the curve between the two converged points. When the curve changes direction (from increasing to decreasing or vice versa) in the interval of interest, it could mean that a discontinuity of s occurs in this region, leading us to divide the corresponding interval in two in order to approximate the values (by using the same values of parameters) in each half. As the approximations are obtained on smaller intervals, the curve is better approximated and it enables us to determine if the change of direction is due to a discontinuity or not. By increasing the number of approximated points, the jump at a discontinuity should become more abrupt, which can be detected by measuring the distance between the points. If the distance of the jump increases significantly, it means that a discontinuity occurs in the interval; if not, it simply indicates that the curve changes direction smoothly. This information, as well as the number of discontinuities in the interval, is then used to determine if a bifurcation point occurs or not in the corresponding interval.
5. Validation of the approach To show that our proposed strategy works well, let us consider again the finite strain incompressible elasticity problem (in the case of initial remeshing), but where this time the results will be compared by using four different approaches. In the first case, the standard approach described in Section 2 will be used. In all other cases, the algorithm of Ligursky and Renard will be incorporated in our procedure, but the parameters and the strategy will differ. For example, in the second case, we will use the same values of parameters as in [20] (dmax ¼ 103 ; dmin ¼ 106 ; dinc ¼ 2; ddec ¼ 0:1). In the third case, we will try increasing these parameters in order to reduce the computational time (dmax ¼ 101 ; dmin ¼ 102 ; dinc ¼ 1:5; ddec ¼ 0:5) while in the fourth case, we will use the same parameters as in the third case, but where further investigation will be activated when needed by looking for changes of direction in the curve of s. The fourth case corresponds to our proposed strategy. Table 1 compares the four approaches by giving the computation time needed to reach a load of F ¼ 600 and the number of bifurcation points detected in all cases. As illustrated in Fig. 3, three bifurcation points should be detected in this time frame. Fig. 11 illustrates the bifurcation diagram in all four cases and shows which bifurcation points were detected. As can be seen, only two approaches led to the detection of all three bifurcation points: our proposed approach as well as Ligursky and Renard’s algorithm with the use of very small parameters. Our proposed approach is however approximately sixty-three times less costly than the algorithm presented in [20], which represents a very significant improvement. We also note that simply using larger values for the parameters in Ligursky and Renard’s algorithm in order to reduce the computational time is not sufficient as only one bifurcation point was detected in that case. Identifying intervals where possible discontinuity of s occur is key to the good performance of the algorithm. Our proposed approach was also tested on a more difficult large deformation elasticity problem. In this example, let us consider the indentation problem introduced in Yamada and Kikuchi [25] where a displacement of magnitude v is applied in the x2 direction on the left half of the top edge of a rectangular geometry of dimension 2 1, as shown in Fig. 12. On the bottom edge, the displacement in the x2 direction is fixed (u2 ¼ 0), while on both side edges, the displacement in the x1 direction is fixed (u1 ¼ 0). For
8
S. Léger et al. / Computers and Structures 191 (2017) 1–11
Fig. 9. Finite strain incompressible elasticity problem: calculated values of obtained using different values for the parameters in Algorithm 1).
s near the first bifurcation point at an approximate load of F ¼ 264 (comparison of the values
Fig. 10. Finite strain incompressible elasticity problem: example of the approximated values of
s when no bifurcations occur in the interval.
9
S. Léger et al. / Computers and Structures 191 (2017) 1–11
Table 1 Finite strain incompressible elasticity problem: computation time and number of bifurcation points detected (comparison of the different approaches for the detection of bifurcation points).
Case Case Case Case
1 2 3 4
(standard approach) (Algorithm 1 with same parameters as in [20]) (Algorithm 1 with larger values for the parameters) (proposed approach)
CPU time (in seconds)
Number of bifurcation points detected
82 70,258 1 103 1109
0 3 1 3
Fig. 11. Finite strain incompressible elasticity problem: bifurcation diagrams (comparison of the different approaches for the detection of bifurcation points).
Fig. 12. Indentation problem [25]: Geometry, boundary conditions and initial mesh.
the simulation, the parameters for the Mooney-Rivlin model will be chosen as k ¼ 40; 000; 000; c1 ¼ 20 and c2 ¼ 0. The initial mesh, which has been remeshed in order to take into consideration the various transitions in the solution from the start, is also shown in Fig. 12. As can be seen, the initial mesh has been refined in the vicinity of point A (which separates two different types of boundary conditions).
By simulating this problem with our proposed approach, a bifurcation point was detected at an approximate load of v ¼ 0:247. By completing the same simulation with the other three approaches, the only other approach capable of detecting this bifurcation point is again Ligursky and Renard’s algorithm with the use of very small parameters (see Table 2). However, the computational cost of our approach is again much smaller. Fig. 13 illus-
10
S. Léger et al. / Computers and Structures 191 (2017) 1–11
Table 2 Indentation problem: computation time and number of bifurcation points detected (comparison of the different approaches for the detection of bifurcation points).
Case Case Case Case
1 2 3 4
(standard approach) (Algorithm 1 with same parameters as in [20]) (Algorithm 1 with larger values for the parameters) (proposed approach)
Fig. 13. Indentation problem: values of
Table 3 Indentation problem: comparison of the six smallest eigenvalues for the bifurcation point detected near v ¼ 0:247. The comparison is made just before (v 1 ) and just after (v 2 ) the detected point. Matrix F 0u
v1
v2
þ0:0009 106
0:0022 106
0:1318 106
0:1305 106
0:1516 106
0:1552 106
0:2574 106
0:2580 106
0:3827 106
0:3827 106
0:3945 106
0:3947 106
CPU time (in seconds)
Number of bifurcation points detected
276 98,236 1 681 2046
0 1 0 1
s near the bifurcation point at v ¼ 0:247.
trates the values of s near this bifurcation point and shows that the detection problems are again due to a discontinuity of s. A discontinuity occurs at a load of v ¼ 0:2415, which is very close to the true bifurcation point. To show that this point is a true bifurcation point, we present in Table 3 the six smallest eigenvalues of the Jacobian matrix F 0u near v ¼ 0:247. We clearly see that the Jacobian matrix has an eigenvalue that passes through zero (which confirms the presence of a bifurcation point). As for the deformed shape related to this computed bifurcation point, it is shown in Fig. 14. We can see that buckling occurs in the vicinity of point A.
The eigenvalue that passes through zero is indicated in bold.
6. Conclusion
Fig. 14. Indentation problem: deformed shape related to the computed bifurcation point (zoom in the singular region).
In this work, we have showed that the detection of bifurcation points is a delicate process and that special attention is needed. When discontinuities of s occur near a real bifurcation point, a standard algorithm for the detection of bifurcation points will most likely not be able to detect the bifurcation point as the change of sign is extremely hard to see. This paper therefore presents a more robust approach for the detection of bifurcation points. It is less costly than the one presented in [20] and enables to detect bifurcation points even when discontinuities of s occur near the bifurcation point. The good performance of the proposed approach was illustrated by numerical examples. This paper also shows that the finite strain incompressible elasticity problem presented in [2,3] exhibits many bifurcation points. As very few large deformation bifurcation problems exist in the literature, the results presented in this work could therefore be used as a benchmark problem to test the efficiency of new algorithms for the detection of bifurcation points.
S. Léger et al. / Computers and Structures 191 (2017) 1–11
Acknowledgements The authors wish to acknowledge the financial support of the Faculty of Science of the Université de Moncton and the Natural Sciences and Engineering Research Council of Canada (NSERC). References [1] Abbot JP. An efficient algorithm for the determination of certain bifurcation points. J Comput Appl Math 1978;4:19–27. [2] Auricchio F, Beirão da Veiga L, Lovadina C, Reali A. A stability study of some mixed finite elements for large deformation elasticity problems. Comput Meth Appl Mech Eng 2005;194(9–11):1075–92. [3] Auricchio F, Beirão da Veiga L, Lovadina C, Reali A, Taylor R, Wriggers P. Approximation of incompressible large deformation elastic problems: some unresolved issues. Comput Mech 2013;52(5):1153–67. [4] Bathe KJ, Dvorkin EN. On the automatic solution of nonlinear finite element equations. Comput Struct 1983;17:871–9. [5] Bathe K-J, Ramm E, Wilson EL. Finite element formulations for large deformation dynamic analysis. Int J Numer Meth Eng 1975;9(2):353–86. [6] Bathe KJ. Finite Element Procedures. New Jersey: Prentice-Hall; 1996. [7] Bois R, Fortin M, Fortin A. A fully optimal anisotropic mesh adaptation method based on a hierarchical error estimator. Comput Meth Appl Mech Eng 2012;209–212:12–27. [8] Bois R, Fortin M, Fortin A, Couët A. High order optimal anisotropic mesh adaptation using hierarchical elements. Euro J Comput Mech/Rev Euro Méc Numér 2012;21(1–2):72–91. 8. [9] Crisfield MA. A fast incremental/iterative solution procedure that handles snap through. Comput Struct 1981;13:55–62. [10] Dhooge A, Govaerts W, Kuznetsov YA. MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs. ACM Trans Math Software 2003;29(2). [11] Riks E. The application of newton’s method to the problem of elastic stability. J Appl Mech 1972;39:1060–6. [12] Georg K. Matrix-free numerical continuation and bifurcation. Numer Funct Anal Optim 2001;22:303–20. [13] Hibbitt HD, Marcal PV, Rice JR. A finite element formulation for problems of large strain and large displacement. Int J Solids Struct 1970;6(8):1069–86.
11
[14] Ikeda K, Yamakawa Y, Tsutsumi S. Simulation and interpretation of diffuse mode bifurcation of elastoplastic solids. J Mech Phys Solids 2003;51 (9):1649–73. [15] Keller HB. Lectures on numerical methods in bifurcation problems. New Delhi: Narosa Publishing House; 1987. [16] Le Tallec P, Vidrascu M. Une méthode numérique pour les problèmes d’équilibre de corps hyperélastiques compressibles en grandes déformations. Numer Math 1984;43:199–224. [17] Léger S, Deteix J, Fortin A. A Moore-Penrose continuation method based on a Schur complement approach for nonlinear finite element bifurcation problems. Comput Struct 2015;152:173–84. [18] Léger S, Fortin A, Tibirna C, Fortin M. An updated Lagrangian method with error estimation and adaptive remeshing for very large deformation elasticity problems. Int J Numer Meth Eng 2014;100(13):1006–30. [19] Léger S, Pepin A. An updated Lagrangian method with error estimation and adaptive remeshing for very large deformation elasticity problems: the threedimensional case. Comput Meth Appl Mech Eng 2016;309:1–18. [20] Ligursky´ T, Renard Y. Bifurcations in piecewise-smooth steady-state problems: abstract study and application to plane contact problems with friction. Comput Mech 2015;56(1):39–62. [21] Neto LB, Machado RD, Hecke MB. Finite element analysis of buried pipelines subjected to buckling. Int J Model Simul Petrol Ind 2009;3(1):23–33. [22] Rheinboldt WC. Numerical analysis of continuation methods for nonlinear structural problems. Comput Struct 1981;13:103–13. [23] Wagner W, Wriggers P. A simple method for the calculation of post-critical branches. Eng Comput 1988;5:103–9. [24] Wriggers P, Simo JC. A general procedure for the direct computation of turning and bifurcation points. Int J Numer Meth Eng 1990;30:155–76. [25] Yamada T, Kikuchi F. An arbitrary Lagrangian-Eulerian finite element method for incompressible hyperelasticity. Comput Meth Appl Mech Eng 1993;102 (2):149–77. [26] Zhou Z, Murray DW. Analysis of postbuckling behavior of line pipe subjected to combined loads. Int J Solids Struct 1995;32(20):3015–36. [27] Zhou Z, Nishida A, Kuwamura H. Applicability of finite element method to collapse analysis of steel connection under compression. Prog Nucl Sci Technol 2011;2:481–5. [28] Zienkiewicz OC, Taylor RL. The finite element method: basic formulations and linear problems, 4th ed., vol. 1. London: McGraw-Hill; 1989.
Contents lists available at ScienceDirect
Computers and Structures journal homepage: www.elsevier.com/locate/compstruc
Improved algorithm for the detection of bifurcation points in nonlinear finite element problems S. Léger ⇑, J. Haché, S. Traoré Département de mathématiques et de statistique, Pavillon Rémi-Rossignol, 18 avenue Antonine-Maillet, Université de Moncton, Moncton E1A 3E9, Canada
a r t i c l e
i n f o
Article history: Received 14 October 2016 Accepted 6 June 2017
Keywords: Finite element method Moore-Penrose continuation method Schur complement Mesh adaptation Algorithm Bifurcations
a b s t r a c t Dealing with bifurcation points when solving large deformation finite element problems is not an easy task. Near such points, the Jacobian matrix becomes singular and the problem becomes difficult to solve numerically. In these situations, increasing heuristically the loading parameter during the simulation in order to follow the solution branch is not an option as this approach usually results in the divergence of the process. Efficient numerical techniques capable of handling the presence of bifurcation points are therefore necessary and continuation methods have proved to be powerful tools when dealing with these kind of issues. In Léger et al. (2015), a new implementation technique based on a Schur complement approach for the Moore-Penrose continuation method, which facilitates the detection of bifurcation points and enables branch following, was presented. This method has proved to perform well in most situations; however, in others (i.e. when mesh adaptation is added to the algorithm), some difficulties appear. In this paper, we therefore present an improved approach, which is much more robust, for the detection of bifurcation points in nonlinear finite element problems. Numerical examples will be presented to show the efficiency of the new approach. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction In recent years, the use of numerical methods to solve industrial problems is constantly increasing. With the complexity of these problems constantly increasing as well, developing robust, efficient and accurate numerical methods is a priority. This is especially true for large deformation hyperelastic problems, where convergence problems are frequently encountered. These problems can be related to different issues, but in most cases they are due to distorted elements that appear in the mesh or to the presence of bifurcation points (where the Jacobian matrix is singular). Bifurcation points, which correspond to eigenvalues of the system, are quantities of interest [26,21,27]. Detecting such points and following the different solution branches passed the critical loads associated with the bifurcation points are thus important and useful to better understand the physical properties of the problem we are solving. The analysis of structural instabilities [6,28], which includes the detection of bifurcation points, is generally based on a numerical continuation method (see [9,4,11,15,23,1,22,24]). In a finite element context, where large systems are frequently encountered, simply looking at the sign of the determinant of the
⇑ Corresponding author. E-mail address: [email protected] (S. Léger). http://dx.doi.org/10.1016/j.compstruc.2017.06.002 0045-7949/Ó 2017 Elsevier Ltd. All rights reserved.
Jacobian matrix for the detection of bifurcation points in not an option as the computation of these determinants end up with machine overflows. A new implementation technique for the Moore-Penrose continuation method, which takes advantage of information already available and facilitates the detection of such points in a finite element context, was therefore introduced in [17]. To validate this new algorithm, the classical elastic beam buckling problem was considered as the bifurcation modes for this problem are well known (see Le Tallec and Vidrascu [16], Ikeda et al. [14]). The finite strain incompressible elasticity problem presented in [2,3] was also used for the validation process. Not only were we able to detect the bifurcation point presented in those papers, but we were also able to detect many more bifurcation points which will be presented in this paper. We note that, at that point, the validation process was done using fairly regular meshes. To solve large deformation problems with the finite element method, two formulations are frequently used in practice: the total Lagrangian formulation [13] which refers to the initial configuration and the updated Lagrangian formulation [5] which refers to the most recently calculated configuration. These two formulations are mathematically equivalent, but the updated Lagrangian formulation tends to perform better in practical applications, as was shown in [18]. It was also shown in [18,19] that remeshing the initial configuration, which does not necessitate any transfer of vari-
2
S. Léger et al. / Computers and Structures 191 (2017) 1–11
ables (meaning that no loss of information results from this), also enhances the performance of the updated Lagrangian method. However, when trying to incorporate the detection of bifurcation points in our complete updated Lagrangian algorithm, we noticed that remeshing had an effect on the detection of bifurcation points. Investigation led us to a better understanding of the phenomena observed and modifications were made to the bifurcation point detection algorithm in order to obtain a much more robust numerical method, which is the object of this paper. The paper is organized as follows. Section 2 describes the implementation technique used for the Moore-Penrose continuation method and explains how this approach facilitates the detection of bifurcation points in a finite element context. Section 3 compares the results obtained when using fairly regular meshes versus adapted meshes. Section 4 describes the improved algorithm for the detection of bifurcation points while Section 5 is devoted to the validation of this improved strategy. 2. Moore-Penrose continuation method When solving large deformation problems, a loading parameter k, corresponding to either external forces or prescribed displacements or both, implicitly drives the deformation. This parameter is typically increased gradually until the desired loading, kmax , is attained. However, increasing heuristically the loading parameter while using a standard Newton method to solve the nonlinear system of equations is highly inefficient as this approach will usually result in the divergence of the algorithm in the neighbourhood of limit points or bifurcation points. A more efficient approach is to use a continuation method where the loading parameter k is explicitly introduced in the system of nonlinear equations, which can be expressed as:
FðxÞ ¼ Fðu; kÞ ¼ 0
ð1Þ
with F : RNþ1 ! RN a smooth function. The vector of unknowns, x, therefore consists of the N degrees of freedom denoted as u plus the loading parameter. The Moore-Penrose continuation method, which was shown to be very robust in the case of large deformation problems (see [18]), is a predictor-corrector method and can be summarized as follows (see [10,17]). Starting from a known point xðiÞ 2 RNþ1 on the solution curve, and given a tangent vector v ðiÞ at that point, the next point xðiþ1Þ 2 RNþ1 on the solution curve as well as the tangent vector at that point v ðiþ1Þ can be obtained using the following algorithm: X 0 ¼ xðiÞ þ hv
AðX
k
Þ dkx
> V k dkx
¼0
AðX k ÞT k ¼ AðX k ÞV k >
Vk Tk ¼ 0 3. Update the solution vector:
X kþ1 ¼ X k dkx 4. Update the tangent vector: k
Z ¼V T
k
6. If kFðX k Þk 6 eF and kX kþ1 X k k 6 ex , convergence attained:
xðiþ1Þ ¼ X kþ1 ;
v ðiþ1Þ ¼ V kþ1
We note that X 0 is the prediction point obtained by making a prediction step of length h in the tangential direction, AðX k Þ ¼ ½F 0u ðX k Þ F 0k ðX k Þ is a rectangular matrix of dimension h i> n ðn þ 1Þ, dkx ¼ dku dkk is a correction vector of dimension ðn þ 1Þ 1; kmax represents the maximum number of iterations allowed while eF and ex are the desired tolerances on F and x respectively. As for T k , it represents a correction on the tangent vector. More details can be found in [17]. As can be seen, to obtain the next point on the solution curve, two systems need to be solved at each iteration. In matrix form, these systems are given by:
F 0u ðX k Þ F 0k ðX k Þ V ku
>
F 0u ðX k Þ F 0k ðX k Þ >
dku
!
¼
dkk
V kk
and
V ku
!
!
T ku
!
T kk
V kk
¼
FðX k Þ 0
!
AðX k ÞV k
ð4Þ
!
0
ð5Þ
where the vectors V k ; T k and dkx were decomposed into their comh i> > ponents in the u and k directions (i.e. V k ¼ V ku V kk ; T k ¼ h i> h i> T ku T kk and dkx ¼ dku dkk ). In a finite element context, constructing and solving these two systems can be very costly. A new algorithm, based on a Schur complement approach, was therefore introduced in [17]. This new approach not only solves the systems by using as much as possible information that is already available, but it can also be incorporated numerically in a finite element code by simply adding a postcondition to the Newton method solver. The following algorithm summarizes the implementation technique of the Moore-Penrose continuation method when using this approach. X 0 ¼ xðiÞ þ hv
ðiÞ
V ¼v For k ¼ 0; 1; 2; . . . ; kmax 0
ðiÞ
(b) Solve for W ku :
ð2Þ
F 0u ðX k Þ W ku ¼ F 0k ðX k Þ (c) Calculate the Schur complement bk :
2. Solve the linear system:
(
Z kZk
F 0u ðX k Þ Dku;k ¼ FðX k Þ
k
¼ FðX Þ
V kþ1 ¼
1. (a) Solve to obtain the standard Newton correction Dku;k (for a fixed value of k):
ðiÞ
V 0 ¼ v ðiÞ For k ¼ 0; 1; 2; . . . ; kmax 1. Solve the linear system:
(
5. Normalization of the tangent vector:
bk ¼ V kk V ku W ku ð3Þ
(d) Calculate the Moore-Penrose correction for the load parameter k:
dkk ¼
1 bk
V ku Dku;k
(e) Apply the Moore-Penrose correction on the standard Newton correction:
3
S. Léger et al. / Computers and Structures 191 (2017) 1–11
dku ¼ Dku;k W ku dkk 2. (a) Calculate the correction for the k component of the tangent vector:
1
T kk ¼
bk
½V ku V ku þ V ku W ku V kk
V ku
¼
þ
W ku
½V kk
3. Numerical results: regular meshes vs adapted meshes Let us consider the finite strain incompressible elasticity problem presented in [2,3], where a square of initial configuration
(b) Calculate the correction for the u component of the tangent vector:
T ku
can be used to follow a specific bifurcation branch (by using it as a prediction direction instead of the usual tangent vector).
T kk
3. Update the solution vector:
X kþ1 ¼ X k dkx 4. Update the tangent vector:
Z ¼ Vk Tk 5. Normalization of the tangent vector:
X0 ¼ ½1; 1 ½1; 1, with coordinates x and y and meshed as in Fig. 1, is considered. For this problem, the top edge of the configuration is traction free, while all other edges are fixed. The body is subjected to a vertical uniform body force F, the incompressible elastic material is described by a Mooney-Rivlin model with parameters k ¼ 40; 000; 000; c1 ¼ 20 and c2 ¼ 0 and P2 P1 elements are used for the computation. By using a standard Newton approach to solve this problem, the algorithm converges to the trivial solution, given by u ¼ 0 and p ¼ Fð1 yÞ, until convergence problems are encountered at a load of approximately F ¼ 265, which corresponds to the critical value detected in [3]. In [17], the algorithm presented in Section 2
Z kZk
V kþ1 ¼
6. If kFðX k Þk 6 eF and kX kþ1 X k k 6 ex , convergence attained:
xðiþ1Þ ¼ X kþ1 ;
v ðiþ1Þ ¼ V kþ1
It was also shown in [17] that this new approach facilitates the detection of bifurcation points. In fact, according to [12], by defining sðxÞ via the system:
JðxÞ B C> 0
YðxÞ sðxÞ
¼
0 ; 1
ð6Þ
where JðxÞ is the Moore–Penrose matrix given by:
JðxÞ ¼
; >
A
v
and where B and C are chosen in such a way that B R ImðJð xÞÞ and C R Im J ð xÞ> (with x the bifurcation point), then sðxÞ will change sign at the bifurcation point x ¼ x. By also solving this additional system using a Schur complement approach, the solution is easily obtained and is given, at each iteration k, by:
Fig. 1. Finite strain incompressible elasticity problem: initial mesh.
bk ¼ V kk V ku W ku (Schur complement) Gkk ¼ b1k ðBk V ku P ku Þ Gku ¼ P ku Gkk W ku h i1 sk ¼ C k Gkk þ C u Gku
0.5
Y kk ¼ sk Gkk
0.3
Gku
0.2 1
0.1
1
where P ku ¼ F 0u ðX k Þ Bu and W ku ¼ F 0u ðX k Þ F 0k ðX k Þ and where the vectors B; C and Y have been decomposed into their components in the u and k directions respectively (i.e. B ¼ ½Bu Bk > , C > ¼ ½C > u C k and Y ¼ ½Y u Y k > ).
As bk and W ku are already known, the detection of bifurcation points can be done by simply solving one more linear system with the matrix F 0u ðX k Þ (to calculate Pku ) and by computing a few scalar products. We note that F 0u ðX k Þ is nothing but the classical Jacobian matrix in a classical Newton method. The detection of bifurcation points therefore only requires adding a few lines of code in the postcondition to the Newton method solver. We also note that while s can be used to monitor bifurcation points, the vector Y
L2
¼ s
k
u
Y ku
0.4
0 -0.1 -0.2 -0.3 -0.4 -0.5 0
200
400
600
800
1000
1200
1400
1600
1800
F Fig. 2. Finite strain incompressible elasticity problem: detected bifurcation points.
4
S. Léger et al. / Computers and Structures 191 (2017) 1–11
Fig. 3. Finite strain incompressible elasticity problem: first six bifurcation modes.
3.5 (a)
3 2.5
(b)
u
L2
2 1.5
(c)
1
(d)
0.5 (e) (f)
0 -0.5 0
200
400
600
800
F
1000
1200
1400
1600
1800
Fig. 4. Finite strain incompressible elasticity problem: bifurcation diagram of the first six bifurcation modes.
was used and a bifurcation point was detected at an approximate value of F ¼ 264, which explains the convergence problems of the standard Newton approach. It was also shown in [17] that another bifurcation point exists for this problem at an approximate value of F ¼ 370. Detecting this bifurcation point was made possi-
Fig. 5. Finite strain incompressible elasticity problem: deformed meshes for the 2nd antisymmetric and symmetric modes of bifurcation.
ble by staying on the fundamental branch passed the first bifurcation point (instead of following the bifurcation branch). At that time, a total body force of F ¼ 500 had been imposed for the computation, leading us to believe that only two bifurcation points existed for this problem. However, by simulating this problem for larger values of F, many more bifurcation points were detected by using our proposed algorithm, as can be seen in Fig. 2. We note that the curve
S. Léger et al. / Computers and Structures 191 (2017) 1–11
Fig. 6. Finite strain incompressible elasticity problem: initial meshes used for the comparison.
that passes through the detected bifurcation points (indicated by the symbol () in the figure) represents the fundamental branch. The bifurcation modes associated with these points, which were validated by looking at the eigenvalues, can also be followed by using the vector Y as described in Section 2. Fig. 3 illustrates the first six bifurcation modes detected for the finite strain incom-
5
pressible elasticity problem and also indicates the approximate load at which these different bifurcations occur. As can be seen, the modes associated with the different bifurcation points alternate between antisymmetric and symmetric. A bifurcation diagram of these six bifurcation modes is shown in Fig. 4, while Fig. 5 illustrates the robustness of our Moore-Penrose continuation method by showing the high level of deformation attained for both the 2nd antisymmetric and symmetric modes of bifurcations. We note that similar results were also obtained for the other modes. Now, as remeshing the initial configuration in an updated Lagrangian context enhances the performance of the numerical method, let us remesh the configuration of Fig. 1 by using the fully optimal mesh adaptation method developed by Bois et al. [7,8] in order to verify that the proposed algorithm still detects the bifurcation points when remeshing is added to the algorithm. We note that for this example, the error indicator not only takes into consideration the displacement. It is also based on the pressure due to the fact that the trivial solution for the displacement is null. By taking into consideration both the displacement and the pressure, we ensure that the error indicator is able to detect the variations in the solution for both eventualities (i.e. when the calculated point lies on the fundamental curve and only the pressure varies as well as when the calculated point lies on a bifurcation branch in
Fig. 7. Finite strain incompressible elasticity problem: calculated values of s near the second bifurcation point at an approximate load of F ¼ 370 (comparison of the values obtained using a regular mesh versus an adapted mesh).
6
S. Léger et al. / Computers and Structures 191 (2017) 1–11
which case the displacement is nonzero). Since the initial mesh is adapted based on the error indicator calculated at the end of the first loading step, we note that both eventualities are feasible depending on the size of the loading step (as well as the step size chosen for the Moore-Penrose continuation method). Adding the remeshing of the initial configuration in our algorithm led us to notice that the detection algorithm did not always work properly. To better understand the reason behind this, let us take a closer look at the values of s calculated during the simulations. For the comparison, we will solve the problem with two different meshes. In one situation, the algorithm will detect the bifurcation point, while in the other, the algorithm will not be able to detect the bifurcation point. Fig. 6 illustrates the two different meshes used for the comparison, while Fig. 7 shows the calculated values of s near the second bifurcation point which is located at an approximate load of F ¼ 370. Before comparing the values of s, let us first remind the reader that while following the solution curve, vector G can become orthogonal to vector C which produces a singularity in the computation of s and indicates that the matrix MðxðsÞÞ, defined by
MðxðsÞÞ ¼
JðxðsÞÞ B C
>
0
ð7Þ
is singular for some value of s. Numerically, the algorithm will see this as a change of sign even though s is not defined there. In [17], we had therefore implemented a second verification test in order to eliminate these false bifurcation points. However, we can easily imagine that difficulties can appear if ever a real bifurcation point is located just after a discontinuity of s (or vice versa). In fact, in that specific case, two changes of sign would happen in a very small interval, meaning that the algorithm could very well not see these two changes of sign, all depending on which points are calculated along the solution curve. By looking at Fig. 7, we can see that the discontinuity of s, which occurred at a load of approximately F ¼ 317 in the case of the regular mesh, now occurs at a load of approximately F ¼ 365 in the case of the adapted mesh, meaning that the discontinuity is now much closer to the real bifurcation point and that the scenario described above could easily happen. To confirm that this is the reason why the bifurcation point has not been detected by the algorithm, we have approximated the values of s at different locations near the bifurcation point. Fig. 8 illustrates these values (indicated by the symbol ()) and confirms our hypothesis. The two consecutive converged points near the bifurcation point are respectively located just before the discontinuity of s and just after the real bifurcation, meaning that both values are positive and that no change of sign is detected. This phenomenon seems to happen more frequently with adapted meshes, or finer meshes, as the discontinuity always seems to take place much closer to the real bifurcation point in those cases. As the algorithm cannot dictate which points are calculated by the continuation method along the solution curve, changes are needed in order to ensure that the bifurcation points are detected no matter which mesh is used initially. This improved and much more robust algorithm is described in the next section. 4. Improved algorithm for the detection of bifurcation points Not many papers deal with the detection of bifurcation points in large deformation problems. Ligursky and Renard [20] have however recently published a paper in which they propose a criterion for detecting bifurcation points as well as a procedure for its realization in a numerical continuation context designed for large problems. Their approach is similar to the one described in Section 2, with the exception that they do not seem to use a Schur
complement approach to solve the different systems. Another difference in their approach is that between each consecutive converged points, they approximate matrix J by using the linear combination:
J :¼ ð1 aÞJðxðiÞ Þ þ aJðxðiþ1Þ Þ where a varies from 0 to 1, xðiÞ is the last converged point and xðiþ1Þ is the current converged point, and then solve the system
J
B
C>
0
Y
s
¼
0 ; 1
ð8Þ
for each a in order to approximate the values of s in between the points. They do not specify that this is to help with the detection of bifurcation points when two changes of sign happen in small intervals, but we can certainly assume that it is, as this procedure gives a better idea of the behaviour of s in between the converged points (as can be seen in Fig. 8 where this procedure was used to approximate the values of s). Their algorithm is described in Algorithm 1. Algorithm 1. Ligursky and Renard’s algorithm [20] for the detection of bifurcation points Require: xðiÞ ; xðiþ1Þ ; v ðiÞ ; v ðiþ1Þ ; dmax > dmin > 0; dinc > 1 > ddec > 0 Step 1: Set nch :¼ 0; a :¼ 0; s0 :¼ 106 ; choose B and C randomly. Step 2: Set
s1 :¼ 106 ; d :¼ dmin , and
J :¼ ð1 aÞJðxðiÞ Þ þ aJðxðiþ1Þ Þ;
M :¼
J C>
B 0
Step 3: Solve
M
Y
s
¼
0 1
Step 4: If ss0 < 0 and js0 j < js1 j, set nch :¼ nch þ 1. Step 5: If js s0 j is large, set d :¼ maxfddec d; dmin g, otherwise if js s0 j is small, set d :¼ minfdinc d; dmax g. Step 6: If a < 1, set a :¼ minfa þ d; 1g; s1 :¼ s0 ; s0 :¼ s and return to Step 2. Step 7: If nch is odd, a bifurcation has been detected. By incorporating this algorithm in our procedure, the bifurcation points are more easily detected. In fact, by using this modified approach (i.e. approximating the values of s in between the converged points), the algorithm now detects the bifurcation point located at F ¼ 370 in the case of the adapted mesh. However, this modification does not completely solve the problem as can be seen in Fig. 9(a) where a bifurcation point is not detected even with the use of the approximated values. In this specific case, a discontinuity of s seems to occur in the interval of interest, but the location of the intermediate points does not enable us to detect the change of sign. By taking smaller values of dmax and dmin , we are able to better approximate the behaviour of s and detect the bifurcation point as can be seen in Fig. 9(b). This however has a fairly high cost associated with it. Instead of taking 0.55 s (as in the first set of parameters) in this specific interval, it takes 3.58 s. As this has to be done between each consecutive converged points, reducing the values of the parameters for the entire simulation would be very costly and
S. Léger et al. / Computers and Structures 191 (2017) 1–11
7
Fig. 8. Finite strain incompressible elasticity problem: approximation of the values of s near the real bifurcation point in the case of the adapted mesh. The approximated values are indicated by the symbol (), while the converged points are indicated by the symbol ().
is not recommended. Furthermore, the number of intermediate points necessary to detect the bifurcation point depends on a number of different aspects (i.e. mesh, problem we are solving, distance between the two converged points, etc.). It is therefore not possible to find optimal values that would work in all cases, unless we take extremely small values for the parameters, which again would result in higher computational times. For example, by using the same values of parameters as in [20], which are extremely small, it takes 514.8 s to detect the bifurcation point in this specific interval. As bifurcations usually only occur a few times during a simulation, there is no need to calculate many intermediate points between each consecutive converged points. However, we do need to calculate a few intermediate points each time as the bifurcation points can easily be hidden when discontinuities of s occur (as was shown in Fig. 8). In all other cases, the intermediate points usually follow a smooth trajectory and lead to the same result as the original algorithm for the detection of bifurcation points, as can be seen in Fig. 10, where in this case, the values of s remain strictly increasing even when more points are calculated. Taking smaller values for the parameters are unnecessary in those cases. Our proposed strategy therefore consists of using larger values of parameters throughout the simulation as to get an idea of the behaviour between the points without generating a high computational cost and then only refining these parameters when needed. To determine if further investigation is needed, we study the shape of the curve between the two converged points. When the curve changes direction (from increasing to decreasing or vice versa) in the interval of interest, it could mean that a discontinuity of s occurs in this region, leading us to divide the corresponding interval in two in order to approximate the values (by using the same values of parameters) in each half. As the approximations are obtained on smaller intervals, the curve is better approximated and it enables us to determine if the change of direction is due to a discontinuity or not. By increasing the number of approximated points, the jump at a discontinuity should become more abrupt, which can be detected by measuring the distance between the points. If the distance of the jump increases significantly, it means that a discontinuity occurs in the interval; if not, it simply indicates that the curve changes direction smoothly. This information, as well as the number of discontinuities in the interval, is then used to determine if a bifurcation point occurs or not in the corresponding interval.
5. Validation of the approach To show that our proposed strategy works well, let us consider again the finite strain incompressible elasticity problem (in the case of initial remeshing), but where this time the results will be compared by using four different approaches. In the first case, the standard approach described in Section 2 will be used. In all other cases, the algorithm of Ligursky and Renard will be incorporated in our procedure, but the parameters and the strategy will differ. For example, in the second case, we will use the same values of parameters as in [20] (dmax ¼ 103 ; dmin ¼ 106 ; dinc ¼ 2; ddec ¼ 0:1). In the third case, we will try increasing these parameters in order to reduce the computational time (dmax ¼ 101 ; dmin ¼ 102 ; dinc ¼ 1:5; ddec ¼ 0:5) while in the fourth case, we will use the same parameters as in the third case, but where further investigation will be activated when needed by looking for changes of direction in the curve of s. The fourth case corresponds to our proposed strategy. Table 1 compares the four approaches by giving the computation time needed to reach a load of F ¼ 600 and the number of bifurcation points detected in all cases. As illustrated in Fig. 3, three bifurcation points should be detected in this time frame. Fig. 11 illustrates the bifurcation diagram in all four cases and shows which bifurcation points were detected. As can be seen, only two approaches led to the detection of all three bifurcation points: our proposed approach as well as Ligursky and Renard’s algorithm with the use of very small parameters. Our proposed approach is however approximately sixty-three times less costly than the algorithm presented in [20], which represents a very significant improvement. We also note that simply using larger values for the parameters in Ligursky and Renard’s algorithm in order to reduce the computational time is not sufficient as only one bifurcation point was detected in that case. Identifying intervals where possible discontinuity of s occur is key to the good performance of the algorithm. Our proposed approach was also tested on a more difficult large deformation elasticity problem. In this example, let us consider the indentation problem introduced in Yamada and Kikuchi [25] where a displacement of magnitude v is applied in the x2 direction on the left half of the top edge of a rectangular geometry of dimension 2 1, as shown in Fig. 12. On the bottom edge, the displacement in the x2 direction is fixed (u2 ¼ 0), while on both side edges, the displacement in the x1 direction is fixed (u1 ¼ 0). For
8
S. Léger et al. / Computers and Structures 191 (2017) 1–11
Fig. 9. Finite strain incompressible elasticity problem: calculated values of obtained using different values for the parameters in Algorithm 1).
s near the first bifurcation point at an approximate load of F ¼ 264 (comparison of the values
Fig. 10. Finite strain incompressible elasticity problem: example of the approximated values of
s when no bifurcations occur in the interval.
9
S. Léger et al. / Computers and Structures 191 (2017) 1–11
Table 1 Finite strain incompressible elasticity problem: computation time and number of bifurcation points detected (comparison of the different approaches for the detection of bifurcation points).
Case Case Case Case
1 2 3 4
(standard approach) (Algorithm 1 with same parameters as in [20]) (Algorithm 1 with larger values for the parameters) (proposed approach)
CPU time (in seconds)
Number of bifurcation points detected
82 70,258 1 103 1109
0 3 1 3
Fig. 11. Finite strain incompressible elasticity problem: bifurcation diagrams (comparison of the different approaches for the detection of bifurcation points).
Fig. 12. Indentation problem [25]: Geometry, boundary conditions and initial mesh.
the simulation, the parameters for the Mooney-Rivlin model will be chosen as k ¼ 40; 000; 000; c1 ¼ 20 and c2 ¼ 0. The initial mesh, which has been remeshed in order to take into consideration the various transitions in the solution from the start, is also shown in Fig. 12. As can be seen, the initial mesh has been refined in the vicinity of point A (which separates two different types of boundary conditions).
By simulating this problem with our proposed approach, a bifurcation point was detected at an approximate load of v ¼ 0:247. By completing the same simulation with the other three approaches, the only other approach capable of detecting this bifurcation point is again Ligursky and Renard’s algorithm with the use of very small parameters (see Table 2). However, the computational cost of our approach is again much smaller. Fig. 13 illus-
10
S. Léger et al. / Computers and Structures 191 (2017) 1–11
Table 2 Indentation problem: computation time and number of bifurcation points detected (comparison of the different approaches for the detection of bifurcation points).
Case Case Case Case
1 2 3 4
(standard approach) (Algorithm 1 with same parameters as in [20]) (Algorithm 1 with larger values for the parameters) (proposed approach)
Fig. 13. Indentation problem: values of
Table 3 Indentation problem: comparison of the six smallest eigenvalues for the bifurcation point detected near v ¼ 0:247. The comparison is made just before (v 1 ) and just after (v 2 ) the detected point. Matrix F 0u
v1
v2
þ0:0009 106
0:0022 106
0:1318 106
0:1305 106
0:1516 106
0:1552 106
0:2574 106
0:2580 106
0:3827 106
0:3827 106
0:3945 106
0:3947 106
CPU time (in seconds)
Number of bifurcation points detected
276 98,236 1 681 2046
0 1 0 1
s near the bifurcation point at v ¼ 0:247.
trates the values of s near this bifurcation point and shows that the detection problems are again due to a discontinuity of s. A discontinuity occurs at a load of v ¼ 0:2415, which is very close to the true bifurcation point. To show that this point is a true bifurcation point, we present in Table 3 the six smallest eigenvalues of the Jacobian matrix F 0u near v ¼ 0:247. We clearly see that the Jacobian matrix has an eigenvalue that passes through zero (which confirms the presence of a bifurcation point). As for the deformed shape related to this computed bifurcation point, it is shown in Fig. 14. We can see that buckling occurs in the vicinity of point A.
The eigenvalue that passes through zero is indicated in bold.
6. Conclusion
Fig. 14. Indentation problem: deformed shape related to the computed bifurcation point (zoom in the singular region).
In this work, we have showed that the detection of bifurcation points is a delicate process and that special attention is needed. When discontinuities of s occur near a real bifurcation point, a standard algorithm for the detection of bifurcation points will most likely not be able to detect the bifurcation point as the change of sign is extremely hard to see. This paper therefore presents a more robust approach for the detection of bifurcation points. It is less costly than the one presented in [20] and enables to detect bifurcation points even when discontinuities of s occur near the bifurcation point. The good performance of the proposed approach was illustrated by numerical examples. This paper also shows that the finite strain incompressible elasticity problem presented in [2,3] exhibits many bifurcation points. As very few large deformation bifurcation problems exist in the literature, the results presented in this work could therefore be used as a benchmark problem to test the efficiency of new algorithms for the detection of bifurcation points.
S. Léger et al. / Computers and Structures 191 (2017) 1–11
Acknowledgements The authors wish to acknowledge the financial support of the Faculty of Science of the Université de Moncton and the Natural Sciences and Engineering Research Council of Canada (NSERC). References [1] Abbot JP. An efficient algorithm for the determination of certain bifurcation points. J Comput Appl Math 1978;4:19–27. [2] Auricchio F, Beirão da Veiga L, Lovadina C, Reali A. A stability study of some mixed finite elements for large deformation elasticity problems. Comput Meth Appl Mech Eng 2005;194(9–11):1075–92. [3] Auricchio F, Beirão da Veiga L, Lovadina C, Reali A, Taylor R, Wriggers P. Approximation of incompressible large deformation elastic problems: some unresolved issues. Comput Mech 2013;52(5):1153–67. [4] Bathe KJ, Dvorkin EN. On the automatic solution of nonlinear finite element equations. Comput Struct 1983;17:871–9. [5] Bathe K-J, Ramm E, Wilson EL. Finite element formulations for large deformation dynamic analysis. Int J Numer Meth Eng 1975;9(2):353–86. [6] Bathe KJ. Finite Element Procedures. New Jersey: Prentice-Hall; 1996. [7] Bois R, Fortin M, Fortin A. A fully optimal anisotropic mesh adaptation method based on a hierarchical error estimator. Comput Meth Appl Mech Eng 2012;209–212:12–27. [8] Bois R, Fortin M, Fortin A, Couët A. High order optimal anisotropic mesh adaptation using hierarchical elements. Euro J Comput Mech/Rev Euro Méc Numér 2012;21(1–2):72–91. 8. [9] Crisfield MA. A fast incremental/iterative solution procedure that handles snap through. Comput Struct 1981;13:55–62. [10] Dhooge A, Govaerts W, Kuznetsov YA. MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs. ACM Trans Math Software 2003;29(2). [11] Riks E. The application of newton’s method to the problem of elastic stability. J Appl Mech 1972;39:1060–6. [12] Georg K. Matrix-free numerical continuation and bifurcation. Numer Funct Anal Optim 2001;22:303–20. [13] Hibbitt HD, Marcal PV, Rice JR. A finite element formulation for problems of large strain and large displacement. Int J Solids Struct 1970;6(8):1069–86.
11
[14] Ikeda K, Yamakawa Y, Tsutsumi S. Simulation and interpretation of diffuse mode bifurcation of elastoplastic solids. J Mech Phys Solids 2003;51 (9):1649–73. [15] Keller HB. Lectures on numerical methods in bifurcation problems. New Delhi: Narosa Publishing House; 1987. [16] Le Tallec P, Vidrascu M. Une méthode numérique pour les problèmes d’équilibre de corps hyperélastiques compressibles en grandes déformations. Numer Math 1984;43:199–224. [17] Léger S, Deteix J, Fortin A. A Moore-Penrose continuation method based on a Schur complement approach for nonlinear finite element bifurcation problems. Comput Struct 2015;152:173–84. [18] Léger S, Fortin A, Tibirna C, Fortin M. An updated Lagrangian method with error estimation and adaptive remeshing for very large deformation elasticity problems. Int J Numer Meth Eng 2014;100(13):1006–30. [19] Léger S, Pepin A. An updated Lagrangian method with error estimation and adaptive remeshing for very large deformation elasticity problems: the threedimensional case. Comput Meth Appl Mech Eng 2016;309:1–18. [20] Ligursky´ T, Renard Y. Bifurcations in piecewise-smooth steady-state problems: abstract study and application to plane contact problems with friction. Comput Mech 2015;56(1):39–62. [21] Neto LB, Machado RD, Hecke MB. Finite element analysis of buried pipelines subjected to buckling. Int J Model Simul Petrol Ind 2009;3(1):23–33. [22] Rheinboldt WC. Numerical analysis of continuation methods for nonlinear structural problems. Comput Struct 1981;13:103–13. [23] Wagner W, Wriggers P. A simple method for the calculation of post-critical branches. Eng Comput 1988;5:103–9. [24] Wriggers P, Simo JC. A general procedure for the direct computation of turning and bifurcation points. Int J Numer Meth Eng 1990;30:155–76. [25] Yamada T, Kikuchi F. An arbitrary Lagrangian-Eulerian finite element method for incompressible hyperelasticity. Comput Meth Appl Mech Eng 1993;102 (2):149–77. [26] Zhou Z, Murray DW. Analysis of postbuckling behavior of line pipe subjected to combined loads. Int J Solids Struct 1995;32(20):3015–36. [27] Zhou Z, Nishida A, Kuwamura H. Applicability of finite element method to collapse analysis of steel connection under compression. Prog Nucl Sci Technol 2011;2:481–5. [28] Zienkiewicz OC, Taylor RL. The finite element method: basic formulations and linear problems, 4th ed., vol. 1. London: McGraw-Hill; 1989.