A neurocomputing model for the elastoplasticity

Comput. Methods Appl. Mech. Engrg. 182 (2000) 177±186

www.elsevier.com/locate/cma

A neurocomputing model for the elastoplasticity q Sun Daoheng a,*, Hu Qiao b, Xu Hao b a

Department of Applied Mechanics and Engineering, Southwest Jiaotong University, Chengdu 610031, People's Republic of China b Mechanical Engineering School of Northeastern University, Shenyang 110006, People's Republic of China Received 21 November 1997; received in revised form 23 July 1998

Abstract To shorten the time of the structural analysis based on the ®nite element method is very important and of interest to engineers and mechanics. In this paper, on the foundation of the ``Parametric Variational Principle for elastoplasticity'', a new variational principle called ``Two-type-variable Variational Principle'' is proposed so that the neural networks theory can be used to solve the elastoplastic problems. For the two-type-variable minimum potential energy principle, the energy function and the constraint conditions of the neural networks are built up, and the structure and parameters are given. As an example, the simulation on the neurocomputation of a simple truss is taken. It is shown that the neurocomputed conclusions are well agreement with that of theoretical analyzed and the solution of the elastoplasticity can be obtained within an elapsed time of the circuit time-constant (nano-second-order). Ó 2000 Published by Elsevier Science S.A. All rights reserved. Keywords: Elastoplasticity; Finite element analysis; Variational principle; Neural networks

1. Introduction To shorten the time of the structural analysis based on the ®nite element method is signi®cance and of interest to mechanics and engineering. Many researches on the problem have been done in last decades. The computation of the elastoplasticity is the foundation of analysis of many complex mechanical behaviors such as fatigue, damage and fracture etc. But, because of the nonlinear relationship between stress and strain in plasticity, whether an element is in elastic or plastic, and in loading or unloading must be assessed according to the loading condition and to the loading criterion, respectively. To solve an elastoplastic problem, it needs iterate many times based on the Deformation Theory of Plasticity (DTP). While on the Flow Theory of Plasticity (FTP), the load have to be divided into many segments and in each of which it needs iterate many times too. The computing process is rather complicate and the time spent is very much [1]. The Parametric Variational Principle (PVP) that is used to solve the mathematical and physical problems with inde®nite boundary was proposed and developed by Prof. W.X. Zhong and coworkers. The relative algorithm of the parametric quadratic programming (QP) was built up [2,3], in which the e€ect of nonlinear on the numerical computation is weaken greatly. Therefore, the computing process is simpli®ed. In Ref. [4], a ®nite element model based on the Neural Networks (NN) theory for elastic problem was proposed and the Modi®ed TH (M-TH) net was given to solve an example. It is shown that the computation can be completed within an elapsed time of the circuit time-constant. In this paper,

q

The subject was supported by the Natural Science Foundation of People's Republic of China. Corresponding author. Present address: Department of Electromechanical Engineering, Xiamen University, Ximen, 361005, People's Republic of China. *

0045-7825/00/$ - see front matter Ó 2000 Published by Elsevier Science S.A. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 9 9 ) 0 0 0 9 1 - 2

178

Sun Daoheng et al. / Comput. Methods Appl. Mech. Engrg. 182 (2000) 177±186

on the foundation of the PVP, a new type of variational principle called Two-type-variable Variational Principle (TVP) is put forward to adapted for neurocomputing the elastoplasticity. The computation of the elastoplasticity is transformed to the QP with inequality constraints, which can be solved by using the neural networks. The solution can be obtained within an elapsed time of the circuit time-constant (ns).

2. Neurocomputation and optimization theory The research on the NN, with the high nonlinear properties and large-scale parallel distributive processing abilities, has made remarkable progress presently. In 1982, the concept of computational energy functional was introduced, the stabilizing criterion was given for the Hop®eld neural networks by John J. Hop®eld [5]. The energy functional's minimum value point of a NN is corresponding to a stable equilibrium state of the system. If a stable attractor of a NN is considered as the local minimum value point of a potential energy function (generalized objective functional), an optimization can be transformed to a stabilizing procedure of a certain NN. With time increasing the potential energy of the NN system decreases until reaching to a stable state. This state is just the minimum value point of the original optimization problem. By using ``simulated annealing method'' or others, the global minimum value point can be obtained. So, the computation may be completed with the system energy ¯owing. In 1986, the ®rst NN used to solve the linear optimum problem was proposed by David W. Tank and John J. Hop®eld, named as TH net. A complicate optimum problem can be solved within an elapsed time of the circuit time-constant in the order of nano second. Hence, a new approach of solving optimum problem was built up [6]. TH net is composed of an objective net and a constraint net, the details can be seen in Ref. [8]. The dynamic equation of ith neuron is m dui ui X ˆ ÿIi ÿ ÿ Aji f …ATj ÿ Jj †; dt Ri jˆ1

…2:1†

Vi ˆ g…ui † i ˆ 1; 2; . . . ; n: The potential energy function of the system can be de®ned as Ep ˆ IV ‡

X X 1 Z Vi F …Aj V ÿ Jj † ‡ gÿ1 …V † dV ; R i 0 j i

…2:2†

which is decreasing with time until the system reaches to stable. So long as the f …x† is selected properly, TH net can be used to solve the constrained optimum problems as follow: min s:t:

I TV Aj V P J :

…2:3†

Modi®ed the connecting mode of the neurons in TH net, it can be used to solve the quadratic programming problem as follow [8,9]: min

1 T V GV ‡ V T I 2

s:t: AV P J :

…2:4† …2:5†

Sun Daoheng et al. / Comput. Methods Appl. Mech. Engrg. 182 (2000) 177±186

179

3. The two-type-variable variational principle for elastoplasticity 3.1. The uni®ed elementary equations There is a structure with domain X. Provided the equilibrium state at some time and the deformation's history are already known, the stress increment and the displacement increment generated next time must satisfy the following equations [4,7]: 1. Equilibrium equation: drij;j ‡ dbi ˆ 0

…in the X†:

…3:1†

2. Strain±displacement relationship: 2deij ˆ dui;j ‡ duj;i

…in the X†:

…3:2†

3. Boundary conditions: drij nj ˆ dpi dui ˆ dui

…on the Sp †;

…3:3†

…on the Su †:

…3:4†

4. Constitutive equations: Stress±strain relation : Loading condition : Flow criterion : where dkp



ˆ0 P0

drij ˆ Dijkl …dekl ÿ depkl †: f …rij ; epij † 6 0:

depij ˆ dkp

…3:5† …3:6†

og ; orij

…3:7†

when f < 0; in elastic or unloading state when f ˆ 0; in plastic loading state

…3:8†

dbi , dpi , dui is the volume force increment, surface force increment and displacement increment on the given boundary, respectively; Dijkl ± elastic coecient tensor; depij ± plastic strain increment; f ± loading function; g ± plastic potential function; dkp ± plastic ¯ow factor. In order to re¯ect directly the e€ect of the loading condition on the stress increment and strain increment, Eq. (3.6) is Tailor expanded with ®rst-order form, and expressed as strain increment [3,4]. At the beginning of a force increment, let f 0 ˆ f …r0ij ; ep0 ij †; then

p 0 f ˆ f …r0ij ‡ drij ; ep0 ij ‡ deij † ˆ f ‡

…3:9† of of drij ‡ p depij : orij oeij

Substituting Eqs. (3.5) and (3.7) to Eq. (3.10) " # of of og of og dekl ‡ dkp ÿ Dijkl : f ˆ f 0 ‡ Dijkl orij oepij orij orij orkl

…3:10†

…3:11†

It means that f is expressed as a linear function of deij and dkp . The yield criterion can be expressed as f …deij ; dkp † 6 0:

…3:12†

Introducing a relaxation variable of m into the inequality (3.12), we can get the equality f …deij ; dkp † ‡ m ˆ 0:

…3:13†

Then Eq. (3.8) can be expressed as dkp m ˆ 0; dkp P 0; m P 0:

…3:14†

180

Sun Daoheng et al. / Comput. Methods Appl. Mech. Engrg. 182 (2000) 177±186

3.2. The two-type-variable minimum potential energy principle (TMPP) In order to apply NN to solve the elastoplastic problems, a new variational principle must be built up. In the PVP, the plastic ¯ow factor dkp is only considered as a controlling parameter rather than variational function. But the dui and dkp can be obtained at the same time by using Lemke or Wolf algorithm. While a NN is applied to solve an optimal problem, its outputs are responding to the independent variables. Therefore, dui and dkp must be considered as independent variables when a new variational principle is built up. We de®ne the potential energy increment functional of a structure X as Z  Pˆ

X

"Z #  Z 1 1 p og Dijkl deij dekl ÿ dk Dijkl dekl dX ÿ dbi dui dX ‡ dpi dui dS : 2 2 orij X Sp

…3:15†

TMPP: In all the possible displacement and the plastic ¯ow factor increment ®elds which satisfy the strain±displacement relationship (3.2) and the displacement boundary condition (3.4), the real solution should be the one which makes the functional P attain to the global minimum value, controlled by the constraint Eqs. (3.13) and (3.14). In the functional P, both dui (or deij ) and dkp are independent. Proof. In TMPP, it is not arbitrary to variate the de and dkp , and it must be constrained with Eq. (3.7). Applied the Lagrange multiplier method to relax the e€ects of the constraint condition, the following equation can be obtained:   Z p p og 0 dX: …3:16† P ˆ P ‡ aij deij ÿ dk orij X The Lagrange multiplier aij can be determined according to the variational fundamental lemma 1 aij ˆ ÿ Dijkl deij : 2

…3:17†

Substituting Eqs. (3.15) and (3.17) into Eq. (3.16), the ®rst variation of P0 will be dP0 ˆ

Z 

 1 1 Dijkl deij d…dekl † ÿ Dijkl depij d…dekl † ÿ Dijkl deij d…depkl † dX 2 2 X "Z # Z dbi d…dui † dX ‡ dpi d…dui † dS : ÿ X

…3:18†

sp

Considering Eq. (3.5), we can obtain "Z # Z  1 drij d…deij † ‡ deij d…drij † dX ÿ dbi d…dui † dX ‡ dpi d…dui † dS dP ˆ X 2 X sp "Z # Z Z 1 d…drij deij † dX ÿ dbi d …dui † dX ‡ dpi d…dui † dS : ˆ X 2 X sp 0

Z

…3:19†

This is just the classical potential variational principle. According to Eq. (3.16), the equation P ˆ P0 : It means that the TMPP is equivalent to the classical potential variational principle.

…3:20† 

Sun Daoheng et al. / Comput. Methods Appl. Mech. Engrg. 182 (2000) 177±186

181

3.3. The two-type-variable minimum complimentary energy principle (TMCP) In TMCP, the complimentary energy increment functional of a system X is de®ned as Z   Z  1 ÿ1 1 og Dijkl drij drkl ‡ dkp drij dX ÿ d ui drij nj dS Pˆ 2 orij X 2 Su

…3:31†

in which, Dÿ1 ijkl is as the inversion of the elastic tensor. TMCP: In all the possible stress and the plastic ¯ow factor increment ®elds which satisfy the equilibrium Eq. (3.1) and the force boundary condition (3.3), the real solution should be the one which makes the functional P attain to the global minimum value, controlled by the constraint Eqs. (3.13) and (3.14). In functional P, both drij and dkp are independent. The proof is omitted. 4. The ®nite element form of the TMPP By using FEM, the continuum X is discretized into many elements, for example, n elements (n1 elastoplastic elements, n2 force boundary elements). Introducing the displacement shape function of an element, we can obtain du ˆ N d;

…4:1†

de ˆ Bd:

…4:2†

In addition, let kp ˆ dkp so that it is expressed simply.Where, d ± displacement increment vector, N ± displacement shape function, B ± nodal displacement±strain matrix …B ˆ L…r†N †. Assembling all the elements and introducing the displacement boundary condition, we can obtain the total potential energy increment of a system and the controlling equations as follow:   1 T T 1 p d Kd ÿ d Uk ‡ q …4:3† min 2 2 8 p < m ˆ ÿQd ‡ U k ‡ d; p T …k † m ˆ 0; : p k P 0; m P 0;

s:t:

…4:4†

where Kˆ

n1 Z X Xe

eˆ1

BT DB dX;

…4:5†

"  T #T n1 Z X og D B dX; Uˆ e or eˆ1 X

qˆ

n Z X eˆ1

e

X

…N T db† dX ‡

n2 Z X eˆ1

" T # n1 Z X of D B dX; Qˆ e or eˆ1 X

Se

…4:6†

 dS; …N T dp†

…4:7†

…4:8†

182

Sun Daoheng et al. / Comput. Methods Appl. Mech. Engrg. 182 (2000) 177±186

"     T  # T n1 Z X of og of og Uˆ D ÿ dX; p e or or oe or eˆ1 X

dˆÿ

n1 Z X Xe

eˆ1

f 0 dX:

…4:9†

…4:10† T

In Eq.(4.4), the constraint condition …kp † m ˆ 0; m P 0; k P 0 embraces the criterion of setting kp when loading or unloading in Eq. (3.8). This is a nonlinear constraint and it is dicult to cope with when we use NN theory to solve the optimization problem. In order to simplify the NN's structure, we try to introduce T the nonlinear constraint …kp † m ˆ 0 into the objective functional. Theorem. Existing a positive M, when l P M, the optimization described with equations   1 T 1 T d Kd ÿ dT Ukp ‡ q ‡ l…kp † m min 2 2  s:t:

m ˆ ÿQd ‡ U kp ‡ d; kp P 0; m P 0

…4:11†

…4:12†

is equivalent to that of Eqs. (4.3) and (4.4). In order to prove the theorem, a lemma is given as follow. Lemma. In the optimization described with Eqs. (4.11) and (4.12), m > 0…m 6ˆ 0†; kp is an optimization variable, there must be a positive number M, when l P M, the minimum value is happened in kp ˆ 0. T

Proof. We only discuss the part of ÿ…1=2†dT Ukp ‡ l…kp † m in Eq. (4.11). Let y1 ˆ l…kp †T m; y2 ˆ ÿ…1=2†dT Ukp , T

T

T

y1 ˆ l…kp † …ÿQd ‡ U kp ‡ d† ˆ l‰…kp † U kp ‡ …kp † …d ÿ Qd†Š;

…4:13†

y1 is a quadratic function of kp and y2 is a linear function of kp . For an engineering problem, in the light of the FTP, the angle of the plastic potential surface gradient and the loading surface gradient is smaller than 90°. If it is the normal ¯ow problem, the two surfaces are coincidence and the angle equals to zero [4]. Hence, the matrix U in Eq. (4.13) is positive, diagonal and symmetric. y1 P 0, because kp P 0; m > 0. According to the characteristics of quadratic function, there must be d ÿ Qd P 0:

…4:14†

Let y ˆ y1 ‡ y2 1 T T T y ˆ l‰…kp † U kp ‡ …kp † …d ÿ Qkp †Š ÿ …kp † UT d: 2

…4:15†

When y attains to the minimum value point, oy 1 ˆ lU kp ‡ l…d ÿ Qd† ÿ UT d ˆ 0; okp 2

…4:16†

 1 kp ˆ ÿ U ÿ1 ‰l…d ÿ Qd† ÿ UT dŠ: l

…4:17†

When d ÿ Qd > 0, there must be a positive number M, make

Sun Daoheng et al. / Comput. Methods Appl. Mech. Engrg. 182 (2000) 177±186

1 M…d ÿ Qd† ÿ UT d P 0; 2

183

…4:18†



kp 6 0, If l P M.  When d ÿ Qd ˆ 0 and UT d 6 0, there is a positive number M satisfying Eq. (4.18) to make kp 6 0, if l P M.  When d ÿ Qd ˆ 0 and UT d > 0, there is a M to make kp 6 0, if M ! 1.  Considering the domain of kp and y being a monotone increasing function when kp P kp , we can  conclude that when kp ˆ 0, the objective function of Eq. (4.11) attains to the minimum value point.  Now, we can prove the theorem. Proof. For arbitrary constant l, the theorem is naturally tenable, if m ˆ 0. If m > 0, it can be seen from the lemma that kp is an optimization variable vector, there must be a positive  number M, when l P M, the optimal solution of the optimization (4.11), (4.12) lies in kp ˆ 0. Therefore, the optimization described with Eqs. (4.11) and (4.12) is equivalent to that of Eqs. (4.3) and (4.4).  5. Neurocomputation of elastoplasticity In order to simplify the NN's structure, introducing the ®rst equation of (4.12) into (4.11), we can get the standard QP: 1 T V GV ‡ V T I 2

min

…5:1†

s:t: AV P J ; where



K Gˆ ÿ…12 UT ‡ lQ† V ˆ fd kp gT ; I ˆ fq

…5:2†  ÿ …12 U ‡ lQT † ; 2lU



ÿQ Aˆ 0

T ÿ ldg ; J ˆ fÿd

 U ; E

T

0g :

The matrix G is positive and symmetric. For the normal ¯ow problem, U ˆ QT .

Fig. 1. The NN used for computation of elastoplasticity.

184

Sun Daoheng et al. / Comput. Methods Appl. Mech. Engrg. 182 (2000) 177±186

The QP described with Eqs. (5.1) and (5.2) can be solved by using M-TH shown in Fig. 1. The M-TH net is composed of the objective net in the left-hand side and the constraint net in the right-hand [1]. In Fig. 1, f … † is the ampli®cation function in the constraint net f …y† ˆ ayl…ÿy†; where



l…ÿy† ˆ

0 y

…5:3†

y P 0; y < 0;

…5:4†

a ± electric resistor, y ± electric current. g… † is the the ampli®cation function of the neuron in the objective net Vi ˆ g…ui † ˆ bui ;

…5:5†

where, b is scalar. The dynamic equations of the NN shown in Fig. 1 are as follow: objective net C

du u ˆ ÿI ÿ ÿ GV ÿ AT u dt R

…5:6†

constraint net u ˆ f …ÿJ ‡ AV †:

…5:7†

The potential function is de®ned as X X1 1 F …ATj V ÿ Jj † ‡ Ep ˆ I T V ‡ V T GV ‡ 2 R j i

Z

Vi

0

gÿ1 …s† ds;

where, F … † is the inde®nite integral of f … † a F …y† ˆ y 2 l…ÿy† 2

…5:8†

…5:9†

 dEp dV  u C dV ˆ I ‡ ‡ GV ‡ AT u ˆ ÿ dt dt R b dt

2

:

…5:10†

Because of C; b > 0 dEp dV dEp 6 0; when and only when ˆ 0: ˆ 0; dt dt dt It is shown from the above equations that the system potential energy is decreasing with time until the system attains to a stable state which is corresponding to a local minimum value of Ep …V †. In Eq. (5.8), the ®rst and second item stand for the objective functional of a quadratic programming, the third for a measurement of disobeying constraint condition, and the forth is an adjustable item. It can be

Fig. 2. A simple truss.

Sun Daoheng et al. / Comput. Methods Appl. Mech. Engrg. 182 (2000) 177±186

185

Table 1 The simulating value and theoretical conclusion p ˆ 3:0  105 N

Theoretical Simulation

p ˆ 4:60305  105 N

p ˆ 5:2  105 N

d1y  E0

k1  E0 k2  E0

k3  E0

d1y  E0

k1  E0 k2  E0

k3  E0 d1y  E0

k1  E0 k2  E0 k3  E0

13048.9 13042.6

0.0 0.0

0.0 0.0

20038.3 20038.2

0.0 0.0

0.0 0.0

0.0 0.0

0.0 0.0

0.38288 0.38096

24633.6 24632.2

46.3361 0.0 46.2843 0.0

proved that the optimal value of the potential functional described in Eq. (5.8) is the solution of the elastoplastic problem described in Eqs. (4.11) and (4.12). 6. The simulation on neurocomputation of a simple truss There is a simple symmetric truss with three elements, subscribed force of P shown as in Fig. 2. The sectional area of each element is F ˆ 1000 mm2 and the length of each elements is l2 ˆ l0 ˆ 100 mm, l1 ˆ l3 ˆ l0 = cos h; h ˆ 30°. The elastic modulus of material E0 ˆ 2:0  105 MPa, yield strength rs ˆ 200 MPa. The yield criterion of each element is de®ned as f …ri † ˆ ri ÿ rs : We can use the NN mentioned above to analyze the elastoplastic behavior. The analysis should be done in neural computer or through dynamic circuit. In this example, we simulated it on the digital computer. In the objective net, there are 4 neurons …d1y ; k1 ; k2 ; k3 †, with ampli®cation gain of b1 ˆ 100:0; b2 ˆ b3 ˆ b4 ˆ 1:0, input resistor Ri ˆ 100 kX, capacitor Ci ˆ 1 pF. In the constraint net, there are 7 neurons with the penal factor a ˆ 5000 X. The values of simulation and theoretical calculation are shown in Table 1. It can be seen from Table 1 that the relative error between the simulation and the theoretical is less than 5/1000. According to the theorem, when a pole is in the critical state of elastoplasticity, ri ˆ rs ; Qdi ÿ d ˆ 0, the coecient of the nonlinear item in objective function should be l ! 1 to obtain the exact solution. However, in the simulation, given l 6 100, the simulation conclusion is well agreement with the theoretical. According to the stabilization time of a dynamic circuit, it will spend time about 10ÿ12 s to analyze the structure shown in Fig. 2 by using neural FEM and it is independent of the structure's complexity. 7. Conclusions (1) Based on the FTP, the FE computation of elastoplastic mechanics was transformed to a quadratic programming problem with inequality constraints, which can be solved by using neural networks. The computation of the problem is equivalent to the dynamic stabilizing procedure of the NN system, the ®nal stable equilibrium point is corresponding to the solution of the mechanical problem. (2) The NN can be mapped into a dynamic circuit. The computation of elastoplasticity can be completed within an elapsed time of circuit time-constant. Therefore, it can be widely used in the ®elds of needing realtime analysis. References [1] J.Q. Jiang, D. Zhang, The Elastoplastic Finite Element Analysis Method (in Chinese), Aeronautical Press, 1990. [2] P. Zen, Probabilistic Fatigue Damage Properties of Materials and Modern Structural Analysis Principle (in Chinese), Science and Technology Literature Press, 1993.

186

Sun Daoheng et al. / Comput. Methods Appl. Mech. Engrg. 182 (2000) 177±186

[3] R.L. Zhang, W.X. Zhong, The numerical solution for PMPEP by parametric quadratic programming, Journal of Computational Structural Mechanics and Applications 4 (1) (1987) 1±11. [4] Sun Daoheng, Hu Qiao, Xu Hao, Neurocomputing theory on ®nite element analysis of solid mechanics (in Chinese), Chinese Journal of Mechanical Engineering 32 (6) (1996) 20±25. [5] J.J. Hop®eld, Neural networks and physical systems with emergent collective computational abilities, PNAS, USA 79(2) (1982) 2554±2558. [6] D.W. Tank, J.J. Hop®eld, Simple neural optimization networks: An A/D converter signal decision circuit and a linear programming circuit, IEEE Trans. CAS 33 (5) (1986) 533±541. [7] R. Wang, Z.H. Xong, W.B. Huang, The Foundation of Plastic Mechanics, Science Press, 1982. [8] L.C. Jiao, Neurocomputation (in Chinese), Xian Electronic Science and Technology University Press, 1993. [9] P.S. Theocaris, P.D. Panagiotopoulos, Neural network for computing in fracture mechanics methods and prospects of applications, Computer Methods in Applied Mechanics and Engineering 106 (1993) 213±228.