Accelerating the nonlinear analysis of hyper-elastic behavior by GMDH-assisted Newton-Raphson scheme

Abstract

This paper presents an accelerated iterative scheme for nonlinear problems. Commonly, analysis of nonlinear behavior is conducted by the Newton-Raphson (NR) method. It is well-known that the number of iterations required depends on the deviation between the “starting point” and the converged solution. In practice, the solution of previous load step is taken as the “starting point”, while the converged solution of the current load step is not known beforehand. Therefore, difficulties or even non-convergence may occur. Recently, it is suggested that a neural network is employed to predict the solution of the current load step. This prediction is then used as the “starting point” for NR scheme. It is expected, that the true converged solution (of current step) is closer to the prediction by neural network than to the solution of previous load step. As a result, the scheme becomes faster due to less iterations. Obviously, any techniques for time-series forecasting can be used. Here, the Group Method of Data Handling (GMDH) is proposed. Loosely speaking, GMDH is a feedforward neural network without backpropagation. Practically, the GMDH-assisted NR scheme should not take longer time than conventional NR scheme. The advantage of GMDH is fast computation; however, the accuracy may be not as high as a network that has backpropagation. Therefore, careful consideration on the construction of GMDH network is needed. In the current work, the performance of GMDH-assisted NR scheme is investigated in analysis of hyper-elastic behavior, which involves both geometrical and material nonlinearity. A study on the influence of activation function on the accuracy is presented. Also, it is found that prediction for incremental displacement (between the current load step and the previous load step) could be better than prediction of displacement of the current load step.

INTRODUCTION

Nonlinear behaviors are usually encountered in engineering problems. Therefore, an efficient algorithm for nonlinear analysis is always needed. In practice, the iterative Newton-Raphson (NR) scheme is the most commonly used, due to simple implementation and quadratic convergence rate. The method has been widely applied in analyses of unsaturated flow 1 , plastic deformation 2 , 3 , geometrical nonlinearity 4 , 5 , 6 , hyper-elastic behavior 7 , 8 , 9 , temperature-dependent heat transfer 10 , 11 , 12 and many other types of problem.

It is common knowledge if the deviation between the converged solution and the “initial guess” or “starting point” is large, difficulties may occur 13 , 14 , 15 , being reflected in the large number of iterations. Non-convergence may even be encountered. Practically, the converged solution of previous load step is used as the “starting point” for the current load step. Usually, small step size is applied, with the hope to increase the possibility of convergence. However, the process would be very time-consuming. In 2001, Kim and Kim introduced the employment of neural network to predict the starting point 13 . A pre-analysis is required to compute parameters that are characteristic to the pattern of the three previous converged solutions. The neural network is then trained to learn the pattern and estimate a starting point for iterations of the current load step. Recently, Nguyen T. N. et al. 14 directly exploits the time-series forecasting of GMDH-type neural network to predict the converged solution of the current load step. The predicted values are then used as the starting point for NR scheme. Ref. 15 further discussed that a careful selection of parameters for GMDH-network could reduce the number of iterations to one iteration. Nevertheless, Refs. 14 , 15 were limited to geometrically nonlinear analysis of shell structures. Furthermore, from a practical point of view, a neural-network-assisted NR scheme should take less time than the conventional NR scheme. Probably this is the reason that there were not much publications on this area, since it is difficult to select a network that satisfies both accuracy and speed.

For analysis of hyper-elastic behaviors of neo-Hookean type materials, which include both geometrical and material nonlinearities, there were some attempts to accelerate the computational process, for e.g. by POD-DEIM 16 , or by the reduced-basis technique namely Combined Approximations (CA) 8 . However, the purpose of Refs. 8 , 16 is saving time in each iteration, rather than reduction of number of iterations as in Refs. 14 , 15 . Being inspired from the Refs. 14 , 15 , in this paper the GMDH-assisted NR scheme is further investigated for analysis of hyper-elastic solids.

This report is organized as follows. Right after the Introduction, a brief review on hyper-elastic behavior is presented in Section 2. The time-series forecasting GMDH network is described in Section 3. Section 4 is reserved for numerical results and discussion. Finally, concluding remarks are drawn in the last section.

HYPER-ELASTIC BEHAVIOR

Let us consider a two-dimensional (2D) solid domain Ω , bounded by Γ , see Figure 1 . It is subject to a body force b , while traction t is applied on the part Γ _t of the boundary, and displacement constraints takes place on Γ _u .

Figure 1 . Sketch of a 2D solid domain

By denoting x the current position of a point in Ω , and X the position of that point in the initial configuration of the body, the displacement vector is u = x – X . Using finite element analysis, the governing equation for equilibrium problem can be written using the initial configuration as follows

where the internal force, F _int , and the external force, F _ext , are given by

Here, S is the second Piola-Kirchhoff stress tensor, N is the vector of shape functions and B ₁ is the derivative operator. Applying the total Lagrangian approach, linearization of Equation (1) for Newton-Raphson iterative scheme takes the following form 8 , 9

where δ u is the change of displacement between two consecutive iterations and K _tan is the tangent stiffness

In Equation (6), D is the fourth-order constitutive tensor, which relates strain and stress components. The matrix B ₁ and B ₂ are calculated by

where n is the number of node within an element. For each local node i ( i = 1, 2, …, n ), we have

In Equation (10), F _IJ are components of the deformation gradient tensor F , which are given by

where 1 is second-order identity tensor. The “comma” sign (,) in Equation (10) and (11) denotes spatial derivative; for e.g., N _i,1 is the derivative of N _i with respect to X ₁ .

The matrix in Equation (7) stores the components of the second Piola-Kirchhoff stress

The second Piola-Kirchhoff stress tensor S (see Equation (2)) is symmetric and can be calculated by

where C = F ^T F is the right Cauchy-strain tensor. I ₁ , I ₂ , I ₃ are the three invariants of tensor C :

The strain energy function W is characteristic to each type of hyper-elastic material. Under the assumption of neo-Hookean materials, W is given by

where κ and μ are the bulk and shear moduli, respectively, and J = det( F ). From Equation (14), the second Piola-Kirchhoff stress tensor S is determined by

The constitutive tensor D is obtained by

METHODOLOGY: GMDH-TYPE NEURAL NETWORK

The group method of data handling (GMDH) 17 is a self-organizing deep learning technique. GMDH is quite similar to a deep neural network, however the number of hidden layers and the number of neurons in each hidden layer is determined by the network itself during the training stage. Another requirement is that output layer of GMDH has only one neuron.

Construction of the first hidden layer

Each neuron of the first hidden layer is equipped by one activation function and one transfer function. The activation function is usually a polynomial that takes k inputs. For example, the activation function in form of a bi-variate ( two inputs ) quadratic polynomials is written as follows

where are the two arbitrary values taken from the input layer. The neuron then produces one output value by the transfer function, which can be chosen as identity function or sigmoid function

Identity function:

Sigmoid function:

Because each neuron takes two inputs, the number of neurons in the first hidden layer can be determined by

where n is the number of inputs.

The neuron is applied to all samples of the training dataset. Assuming that there are s samples in the training dataset , the following equations are obtained for each neuron

Or in the matrix form

where X ( size s -by-6 ) contains the input values and A (size 6-by-1) stores the coefficients. Requiring that produced values from transfer function are equal to the true output of the dataset, we have

The vector of coefficients, A , is then calculated by

Once all the neurons of the layer have been constructed, the mean squared error (MSE) between the values produced by each neuron with the targeted output (from the validation dataset ) can be evaluated. The neurons are then sorted by MSE (in ascending order). A threshold can be applied here to remove neurons that have MSE higher than threshold.

Construction of other hidden layers

The construction of other hidden layers is similar to that of the first layer, except that the inputs for the r ^th hidden layer are taken from the produced values from the previous layer, i.e. the ( r -1) ^th hidden layer.

The construction stops if: (i) the maximum number of hidden layers (specified by the user) is reached, or (ii) the MSE of the best neuron of the current hidden layer is not lower than the MSE of the best neuron of the previous layer. Case (ii) indicates that the result cannot be improved. Therefore, the current layer will be removed and the construction process is terminated.

The final output of the GMDH network is the output value of the best neuron of the very last layer.

GMDH-assisted Newton-Raphson scheme

Figure 2 presents an illustration for Newton-Rapshson (NR) algorithm. Typically, in the beginning of the current load step, namely load step ( t +1), the converged solution of the previous step, i.e. load step ( t ), is taken as the “starting point” for iterations (the solid dot in Figure 2 ). Here, following Refs. 14 , 15 , the GMDH network for time-series forecasting is employed to predict the converged solution of the current step. That predicted value is then used as the “starting point” for NR scheme. It is expected that the predicted value will be closer to the true converged solution and thus, reduction in the number of iterations can be achieved.

Assuming that M solutions from load step ( t - M ) to load step ( t ) have been known. Given the number of delays (i.e. the number of inputs in each sample), the data can be arranged into input sets and target sets. Figure 3 is an illustration of data preparation for M = 8 and delays = 2, where there are 7 samples of data (each sample has 2 inputs and 1 targeted output). These data are sent to GMDH to train and predict the solution of load step ( t +1). A portion of samples are used for training (training dataset) and the rest are for validation (validation dataset).

Figure 2 . Illustration of Newton-Raphson scheme

Figure 3 . Illustration of data preparation for GMDH network

It is noted that similar to any other deep learning technique, normalization of data could be necessary. Furthermore, a GMDH network is needed for each unknown degrees of freedom.

In this paper, the performance of GMDH-assisted NR scheme for analysis of hyper-elastic behavior is investigated. Unlike in previous works 14 , 15 , where only geometrical nonlinearity is considered, both geometrical and material nonlinearities are involved in behavior of hyper-elastic solids. Therefore, the problem is more complicated.

It is noted that the components of tensor D (see Eq. (6) and Eq. (20)) at an arbitrary point are dependent on the strain values at that point, indicating material nonlinearity. In Refs. 14 , 15 , where material nonlinearity is not considered, tensor D would be constant. Here, there are two sources of nonlinearities, meaning that complexity is increased. As a result, it is more difficult for the neural network to learn and provide prediction with high accuracy.

RESULTS AND DISCUSSION

The proposed model is applied to study behavior of a curved beam, as sketched in Figure 4 . The beam is subject to an inclined uniform load (45 ^o ) at one end, while the other end is fixed. The material is assumed to be neo-Hookean type with bulk modulus κ = 120.291 MPa and shear modulus μ = 80.194 MPa. The uniform load q = 0.5 N/mm ² . The problem domain is uniformly discretized by 9 x 100 8-node quadrilateral elements (9 elements along the radial direction and 100 elements along the circumferential direction). There are 2919 nodes in total.

Figure 4 . Sketch of the curved beam problem

Figure 5 . Number of iterations required by Conventional NR for the curved beam problem

Figure 6 . Number of iterations required by GMDH-assisted NR for the curved beam problem

The load is gradually increased by 40 steps, i.e. incremental size Δq = 0.0125 N/mm ² . The conventional NR scheme is conducted to solve for the first M = 9 load steps. The rest of load steps are aided by GMDH to estimate the starting point. By default, identity function is used as transfer function and number of delays is 3.

Comparison between conventional NR and GMDH-assisted NR

Figure 5 and Figure 6 respectively present the number of iterations in each load step, for conventional NR and GMDH-assisted NR using tri-variate, quadratic polynomials (in short “3-quadratic”) as activation function. It is observed in Figure 6 that from load step 10 to load step 40, with the aid of GMDH, the number of iterations in each load step is generally reduced, resulting a total of 141 iterations, which is much less than 207 iterations needed by the conventional NR.

Table 1 Comparison of elapsed time between Conventional NR and GMDH-assisted NR

Table 2 Comparison of elapsed time variations of GMDH-assisted NR

Figure 7 . The curve of vertical reaction force and vertical displacement at point A

The reduction in number of iteration indeed boosts the computational speed, as elapsed time in GMDH-assisted NR (The time for running GMDH is already included) is much less than that in conventional NR, see Table 1 . Roughly 30% of total time can be saved. In Ref. 15 , even one iteration for each load step can be achieved, when only large deformation (geometrical nonlinearity) is considered. The hyper-elastic behavior is more complicated in nature, since the nonlinearities come from both large deformation and the stress-strain relation.

Figure 7 depicts the curves of vertical reaction force and vertical displacement at point A (see Figure 4 ). The results by GMDH-assisted NR and Conventional NR is in good agreement. This is as expected, because the overall procedure of NR is not altered. The role of GMDH is simply providing a better starting point for iterations.

Effects of activation function

Next, the effect of activation function is investigated. Four types of polynomials are taken into account: bi-variate quadratic (2-quadratic), tri-variate quadratic (3-quadratic), bi-variate cubic (2-cubic) and tri-variate cubic (3-cubic). Among the four variations, 3-quadratic is computationally the most efficient, as reported in Table 2 . On the other hand, non-convergence occurs in 2-cubic.

Prediction of incremental displacement

Here, we predict the incremental displacement between load step ( t ) and load step ( t +1), Δ u = u _t+1 - u _t , instead of directly predict the displacement at load step ( t +1), u _t+1 . The same four variations as in sub-section 4.2 are considered.

The performance of GMDH in Table 3 is in general better than that in Table 2 , in terms of number of iterations. Convergence is achieved by all four variations. In both Table 2 and Table 3 , the tri-variate polynomials are better choice for activation function than the two-variate polynomials.

Table 3 Comparison of elapsed time variations of GMDH-assisted NR, in which the incremental displacements are predicted

CONCLUSIONS

The GMDH-assisted NR has been successfully further extended for analysis of hyper-elastic behavior of two-dimensional solids. The prediction by GMDH provides a better starting point for NR iterative scheme, such that the number of iterations can be reduced. As a result, a large amount of computational time is saved.

The efficiency of the proposed scheme comes from the quick process of GMDH. This is important, because the training is online, i.e. it is conducted when the problem is solved. There is no pre-training. Furthermore, several GMDH networks are needed every load steps. The number of GMDH networks is equal to that of unknown degrees of freedom. Therefore, the necessary to have fast computation in each individual network is more pronounced.

It is aware that the accuracy of prediction by GMDH is crucial. The choice of transfer function and activation function would have influence on the accuracy. A comparative study on activation function, while identity function is selected as transfer function, has shown that tri-variate polynomials would result in better performance than two-variate polynomials. Currently, only polynomials are considered for activation function. The possibility of other types of activation function should also be studied in future works. Even the GMDH could be replaced by any other time-series forecasting network. From practical point of view, an accelerated NR scheme is only beneficial if it is faster than conventional NR scheme. Therefore, any attempts to improve accuracy of prediction should always pay attention to the elapsed time. Furthermore, it is found that prediction for incremental displacement could be possibly more robust than direct prediction for the converged value of displacement.

Last but not least, by a good estimation of starting point, GMDH could help to reduce the number of iterations, but it does not have any role in the computational time of each iteration. On the other hand, the reduced basis approach 8 is efficient to accelerate each iteration but cannot reduce the number of iterations. Therefore, combination of two techniques is promising for future works.

List of abbreviations

GMDH: Group method of data handling

NR: Newton-Raphson

MSE: mean squared error

Conflict of interest

There is no conflict of interest.

Contribution of each author

Minh N. Nguyen contributes every aspect of the manuscript: ideas brainstorming, data checking, writing, proofreading and editing of the manuscript.