行政院國家科學委員會專題研究計畫 成果報告
降伏厚面開關機制之理論實驗與應用(2/2)
計畫類別： 個別型計畫 計畫編號： NSC93-2211-E-002-016- 執行期間： 93 年 08 月 01 日至 94 年 07 月 31 日 執行單位： 國立臺灣大學土木工程學系暨研究所 計畫主持人： 洪宏基 報告類型： 完整報告 處理方式： 本計畫可公開查詢中 華 民 國 95 年 3 月 13 日
行政院國家科學委員會專題研究計畫報告
降
伏厚面開關機制之理論實驗與應用
(2/2)Switching mechanism of yield thick-surface: theory, experiment and applications NSC 93-2211-E-002-016 2004 年 8 月 1 日至 2005 年 7 月 31 日 洪宏基 教授 [email protected] 國立台灣大學土木工程學系 摘要 本報告探討開關機制之理論, 以描述彈塑性及摩擦之實驗現象, 目的在於應用於彈塑性組成建模中的彈、 塑兩相切換, 以及摩擦振動控制的定著、 滑動兩相切換。 研究對象以降伏厚面開關機制為主, 惟仍探討其他可 能性, 如多套互補三元 (多個開關、 塑性當量、 內在時間尺度)、 橫貫臨界分叉等。 首先我們由各種面向研究互 補三元, 探討各種相關提法。 為解決反向載重時應力應變曲線過方, 以及循環載重時回應力循環振盪等傳統彈塑性組成律的缺陷, 本 報告推廣互補三元, 得到多互補三元、 連續互補三元。 互補三元是降伏面開關機制, 而多互補三元是多降伏面 開關機制, 連續互補三元則為連續無窮多降伏面開關機制。 同樣為解決上述缺陷, 本報告又提出了另一種新的 表示法, 將降伏面由零厚度增厚為有厚度, 以增大塑性相的組成方程式的適用空間。 為彌補新表示法因增厚而 喪失的一致性條件, 降伏厚面必須有異於傳統降伏面的開關機制, 本報告探討其公設體系。 此外, 本報告也推廣控制動態系統理論的狀態方程式、 輸出方程式, 得到非常一般的彈塑性演化及狀態方 程式, 包含彈性關係式、 彈塑分解耦合、 塑流及硬軟化規則等。 此非常一般的彈塑性演化及狀態方程式, 配合降 伏面開關機制或降伏厚面開關機制, 可將彈塑性組成理論置於控制動態系統理論之上, 繼續發展、 深化, 並建 立層級組織架構。 關鍵詞: 開關機制、 兩相切換、 互補三元、 降伏厚面、 橫貫臨界分叉、 彈塑性、 摩擦。
Switching mechanism of yield thick-surface: theory, experiment and applications
Hong-Ki Hong
Department of Civil Engineering National Taiwan University
Taipei, Taiwan Abstract
The project was proposed to study mechenisms of switching, with the aim of applying to switching between plastic and elastic phases in the problem of elastoplastic constitutive modeling and switching between sliding and sticking in the problem of vibration control by friction. The focus is on the switching mechanism associated with a yield thick-surface, but other feasible approaches such as multiple complementary trios and transcritical bifurcation are also explored. In particular, a variety of aspects of so-called complementary trios and variants are investigated.
In order to deal with the problems that the models with a conventional yield surface tend to predict a loop of stress-strain curve which is over square in reverse loading and whose back stress is over oscillatory in cyclic loading, we generalize a complementary trio to multiple complementary trios, furthermore to a continuum of complementary trios. When associated with a yield surface a complementary trio acts as a switching mechanism of single yield surface. By analogous association, the multiple complementary trios thus act as a multiple yield surface switching mechanism, and the complementary trio continuum a distributed yield surface switching mechanism. Alternatively, we may thicken the yield surface so that the plastic-phase constitutive equations can develop their capabilities in more ample space in stress space. To remedy the loss of consistency condition due to the thickening, the yield thick-surface must be associated with a switching mechanism different from that associated with the conventional zero-thickness yield surface.
Referring to the state and output equations in control dynamical system theory, we obtain a general system of equations for elastoplatic evolution and state, which embrace most ingredients of elastoplasticity such as elastic relations, elastic-plastic decomposition, plastic flow, hardening-softening, etc. The significance of doing so is that the constitutive theory of elastoplasticity thus has a clear relation with the theories of control dynamic systems and nonlinear ordinary differential equations. On the basis of the latter we can deepen our research on plasticity.
Keywords: switching mechanism, two alternative phases, complementary trio, yield thick-surface, transcritical bifurcation, elastoplasticity, friction.
1 Introduction
For reverse loading and yielding the stress-strain curves predicted by the models with a con-ventional yield surface seem to be over “square” at the elastic-plastic transition. In fact, the stress-strain curves obtained in tests of many metallic materials look very smooth not only in the vicinity of the so-called “yield point” but actually very smooth almost over the whole decending curve (including unloading, reverse loading and yielding down to another unloading). In other words, the reality is that the smoothness is a global phenomenon rather than a local one. For cyclic plasticity the back stress (that is the center of the yield surface) and the yield surface pre-dicted by the models with a conventional yield surface seem to be over oscillatory.
To explore the source of the over-square, let us examine a simple model — the Prandtl-Reuss equation: 1 2G ds dλ + 1 2τy s = de dλ (1)
and compare it with the Maxwell model of viscoelasticity: 1 2G ds dt + 1 2τy s = de dt, (2)
where e and s are the deviatoric tensor of strain and stress, respectively, λ is a measure of plasticity, t is time, and G and τy are material constants. It can be observed that the two equations are almost
identical except that the role of t in viscoelasticity is played by λ in elastoplasticity. We know that when excited cyclically the Maxwell model of viscoelasticity exhibits a very smooth stress-strain loop, but the Prandtl-Reuss model of elastoplasticity instead generates an extremely square loop. Why do the two models governed by analogous constitutive equations generate so different responses? To answer it, let us first take a look at their solutions. The solution to equation (1) is
s(λ) = exp−G(λ − λi) τy s(λi) + 2G Z λ λi exp−G(λ − λ 0) τy de(λ0), (3)
and that to equation (2) is s(t) = exp−G(t − ti) τy s(ti) + 2G Z t ti exp−G(t − t 0) τy de(t0), (4)
where the subscripti and the prime0 are used to indicate the initial value and the dummy variable
he may tend to conclude that both the response curves must be smooth. In fact, he is correct but the premise is that the initial value must be something different from the ultimate value of the exponental process involved. It is the difference (or room or space) which is essential for the exponential process to manifest itself. Indeed, the lack of the difference (or room or space) prevents the Prandtl-Reuss equation to generate smooth response curves although it holds the capability. On the contrary the Maxwell equation has ample space to show off its capability so as to exhibit smooth response curves.
To deepen our understanding of this point, let us give the Prandtl-Reuss equation a mechnical-element explanation. It is known that the Maxwell model of viscoelasticity can be viewed as a series of an elastic spring and a viscous damper where a piston-cylinder construction has been used to represent the viscous damper. The constitutive law governing the viscous damper is
s = 2τy
dev
dt . (5)
In contrast with this t-rate damper law the counterpart for the Prandtl-Reuss equation is a λ-rate damper law
s = 2τy
dep
dλ. (6)
in case dλ > 0. However, the Prandtl-Reuss equation (1) is by no means the whole story of the Prandtl-Reuss model of elastoplasticity; it must be accompanied by the complementary trio:
r 1 2s · s ≤ τy, ˙λ ≥ 0, ˙λ r 1 2s · s = τy˙λ. (7)
Here a superimposed dot denotes differentiation with respect to time, d/dt, and a dot between 2 tensors of the same order in the same dimensional space indicates taking their inner product. In other words, in addition to the piston-cylinder construction for simulating equation (6), the Prandtl-Reuss “damper” is also equipped with an on-off switch, which may be sketched as a (Saint Venant) “friction-slide” for simulating the trio (7). Moreover, even the piston-cylinder construction herein is not a usual one, it is governed by the λ-rate law (6) which already presumes dλ > 0. The complementary trio forces the Prandtl-Reuss equation (1) to be possibly valid only on the yield surface
q
1
2s · s = τy, the thickness of which is zero and which is therefore of zero volume in
stress space. As a consequence, there are little room or space for the exponential process to occur so that all the Prandtl-Reuss “damper” can do is to “slide” without the exponential “relaxation”. On the other hand, the Maxwell equation of viscoelasticity is valid in the whole stress space and there are ample room or space for the exponential process to develop.
To solve the problem raised thoroughly we are thus led to modify the concept of the yield surface to a great extent, allowing it to have thickness such that the elasoplastic equation (e.g., the Prandtl-Reuss equation) has a suitable space (or thickness or volume) to develop its capability of (exponential) “relaxation” in addition to the “friction-sliding”.
2 Complementary trio
We call the following expressions
a ≥ 0, b ≥ 0, ab = 0, (8)
the complemantary trio. It can be interpreted as b = 0 if a > 0,
b ≥ 0 if a = 0,
a is not allowed to be < 0.
(9)
This interpretation distinguishes two branches in the ab plane: the positive a-axis and the nonneg-ative b-axis. In the positive a-axis a > 0 and b = 0, whereas in the nonnegnonneg-ative b-axis b ≥ 0 and a = 0. These features enable the complementary trio to play the role of a swithching mechanism. (The roles of a and b can be reversed and there thus results in another interpretation.)
For elastoplasticity the positive a-axis may be identified with the plastic phase and the non-negative b-axis the elastic phase. For friction the positive a-axis may be identified with the sliding phase and the nonnegative b-axis the sticking phase.
3 Switching mechanism for friction
Denote the interface friction force by r. The following expression [3]
r = (
ry if ˙x > 0,
−ry if ˙x < 0
(10)
is usually employed to represent the two-valuedness of Coulomb’s perfect dry contact friction, where ry is a positive constant to be determined experimentally. If the contact surfaces are horizontal
on the earth ground, ry = µmg, where g is the acceleration of gravity and µ is the coefficient of
friction. This formalism is correct but imcomplete. In fact, the friction force r may take any value between ry and −ry when the interface relative velocity ˙x = 0; therefore, the following expression
provides a more precise description, r = ry if ˙x > 0, ∈ [−ry, ry] if ˙x = 0, = −ry if ˙x < 0. (11)
Notice the distinction between the two-valuedness of r in equation (10) and the multiple-valuedness of r in equation (11). It can be seen that there are three phases in equation (11).
Nevertheless, the formalism (11), although correct, is not complete yet, since it is still lack of a two-way relation between ˙x and r. For completeness we need a flow rule (12) and the complementary trio (13)-(15) as follows:
˙x = r ˙ Λ r2 y , (12) | r |≤ ry, (13) ˙ Λ ≥ 0, (14) | r | ˙Λ = ryΛ,˙ (15)
where ˙Λ is the friction power, so that Λ is the dissipated energy due to friction. It is not difficult to verify that equations (12)-(15) imply equation (11), but conversely that equation (11) does not suffice to lead to equations (12)-(15). The proof can be found in Ref. [3].
A significant merit of the complementary trio can be observed by comparing equations (12)-(15) which take advantage of the complementary trio with equation (11) which does not take advantage of the complementary trio. It is difficult to extend the one-dimensional model of equation (11) to higher dimensions. However, generalization of the model (12)-(15) to higher dimensions is rather straightforward as shown below:
˙x = r ˙ Λ r2 y , (16) krk ≤ ry, (17) ˙ Λ ≥ 0, (18) krk ˙Λ = ryΛ,˙ (19)
where krk is the Euclidean norm of r and r and x are the n × 1 matrices of friction force and displacement, respectively. It is clear that equation (16) may be replaced by
where ˙λ is a proportional multiplier with λ = ` ry = Λ r2 y , (21)
` being the sliding length and Λ the dissipated energy due to friction. From equations (18) and ry > 0 we have ˙λ and ˙` ≥ 0, which means the sliding length is never decreasing, an obvious yet
essential fact.
It is seen that the complementary trio (17)-(19) plays the role of switching nechanism for Coulomb’s dry friction.
4 Switching mechanism for elastolasticity
Let us consider a constitutive model of perfect elastoplasticity as postulated below [4].
˙q = ˙qe+ ˙qp, (22) ˙ Q = k ˙qe, (23) Q ˙λ = Qy˙qp, (24) kQk ≤ Qy, (25) ˙λ ≥ 0, (26) kQk ˙λ = Qy˙λ. (27)
Here Q = (Q1, Q2, · · · , Qn) denotes the n-dimensional generalized stress vector and ˙q = ( ˙q1, ˙q2,
· · · , ˙qn) denotes the corresponding generalized strain rate vector. The above constitutive model is
re-postulated from the celebrated Prandtl-Reuss equation formulated by Prandtl (1924) and Reuss (1930), which is well known as the simplest (yet rather useful) three-dimensional constitutive law for describing a class of linearly elastic-perfectly plastic materials.
It is seen that the complementary trio (25)-(27) plays the role of switching nechanism for perfect elastoplasticity. For (more general) plasticity the complementary trio (8) reads
f ≤ 0, ˙λ ≥ 0, f ˙λ = 0, (28)
where f = 0 is the yield surface and λ is the measure of plasticity. The complementary trio (28) offer a switching mechanism for the models of (more general) plasticity.
Let a, b in the complementary trio a ≥ 0, b ≥ 0, ab = 0 be dimensionless, and define c := a − b√ 2 , d := a + b √ 2 ,
which are obviously dimensionless. Thus, the complementary trio a ≥ 0, b ≥ 0, ab = 0 can be simply expressed by
d =| c |, or
a + b =| a − b | .
We may further write it in the form φθ(a, b) := θ(| a − b |) − θ(a) − θ(b) = 0,
in which θ(a) = a.
However, the θ function may be generalized to be any strictly increasing function [7] θ : R → R with θ(0) = 0.
The φθ function is one of the so-called NCP-functions which serve to transform a complementarity
or related problem into a nonlinear and frequently nondifferential system of equations or into a minimization problem [2].
The a, b ∈ R in the complementary trio are real numbers. They may be generalized to N -dimensional real vectors a = (a1, · · · , aN), b = (b1, · · · , bN) ∈ RN and the N -fold complementary
trios:
aν ≥ 0, bν ≥ 0, aνbν = 0 ∀ν = 1, · · · , N, (29)
or, written in vector form,
a ≥ 0, b ≥ 0, a · b = 0. (30)
As usual, the symbols ≥ etc. in vectorial inequalities are understood to mean so componentwise. The equivalence of the multiple complementary trios (29) and the vector complementary trio (30)
is established as follows: (29)3 to (30)3 is obvious by addition, while (30) to (29) may be proved
by contradiction. Let I := {1, · · · , N } = I6=∪ I= and I6=∩ I= = ∅. Suppose
aνbν 6= 0 ∀ν ∈ I6=, (31)
aνbν = 0 ∀ν ∈ I=. (32)
From (30)1 and (30)2, aνbν ≥ 0 ∀ν ∈ I. This and (31) ensure aνbν > 0 ∀ν ∈ I6=. Combining this
with (32) gives a · b > 0, which contradicts (30)3 unless I6= = ∅ and I= = I. This completes the
proof.
To describe the multiple complementary trios (29) it is often said that at least one variable vanishes in each pair of corresponding sign-constrained variables. To describe the vector comple-mentary trio (30) it is often said that the vectors a and b are sign-constrained (or sense-constrained) and orthogonal.
The N -fold complementary trios may be furthermore made to be continuous as follows: a(ν) ≥ 0, b(ν) ≥ 0, a(ν)b(ν) = 0 ∀ν ∈ I ⊆ R, (33) or a(ν) ≥ 0, b(ν) ≥ 0, Z I a(ν)b(ν)dν = 0 ∀ν ∈ I ⊆ R.
6 Multiple switching mechanism
Mimicking the N -fold complementary trios (29), we may devise the N -fold switching mecha-nism as follows:
fν ≤ 0, ˙λν ≥ 0, fν˙λν = 0 ∀ν = 1, · · · , N. (34)
The following is a model with yield polyhedron and interacting linear hardening: ˙q = ˙qe+ ˙qp, (35) Q = Qa+ Qb, (36) ˙qe = C ˙Q, (37) ˙qp = N ˙λ, (38) NTQ˙b = H ˙λ, (39) NTQa≤ Qy, (40) ˙ λ ≥ 0, (41) QTaN ˙λ = QTyλ,˙ (42)
where T denotes matrix transpose. In the model there are four constant matrices of material properties. C is the elastic compliance matrix of order n × n. N is the yield polyhedral normal matrix of order n × N , each column being a unit vector normal to one side of the yield polyhedron. Each side is a yield hyperplane representing a yield mode. H is the kinematic modulus matrix of order N × N , whose off-diaginal entries indicate hardening interactions between different yield modes. Qy is the initial yield stress matrix of order N × 1, each entry representing the initial yield
stress of a yield mode.
The formulation of the preceding model may be rearranged as follows:
˙q = C ˙Q + N ˙λ, (43)
7 A model with two complememtary trios and two measures of plasticity
A two-complementary-trios cyclic elastoplasticity model may be postulated as follows: ˙eij = ˙eeij + ˙e
p ij, (45) sij = sdij + scij + saij, (46) ˙sij = 2G ˙eeij, (47) ˙epij = saij3 ˙λa 2ha + scij3 ˙λc 2hc , (48) ˙saij = −C1˙λcsaij − 2 3k 0 a+ γcC2 ˙epij − C3˙λcepij + (1 − γc) ˙sij+ C4˙λcsij, (49) ˙scij = −C1˙λcscij + 2 3k 0 c+ γcC2 ˙epij + γc˙sij, (50) fa ≤ 0, (51) ˙λa ≥ 0, (52) fa˙λa = 0, (53) fc ≤ 0, (54) ˙λc ≥ 0, (55) fc˙λc = 0, (56) in which eij, eeij, e p ij, s a
ij, scij, sdij and sij are respectively the deviatoric tensors of strain, elastic
strain, plastic strain, active stress, contact stress, eccentric stress and stress, fa:= r 3 2s a ijsaij − ha (57)
is the von Mises yield function, and
fc := r 3 2s c ijscij − hc (58) is a contact function.
The material functions are dependent on the equivalent plastic strain ¯ep defined by
¯ ep = Z t t0 ˙¯ ep(ξ)dξ, e˙¯p := 2 3˙e p ij˙e p ij 12 ≥ 0. (59)
It is seen that ¯ep is a time-like parameter because its time derivative ˙¯ep is nonnegative, and hence is an indicator for the irreversible change of material properties.
The above ˙λc is defined as ˙λc := γce˙¯p (60) with γc := sc ij q 2 3hc sa ij q 2 3ha −h 0 c k0 c . (61)
The material functions h0a, ka0, h0c and k0c are assumed to be related by
(h0a+ ka0) 1 − exp " a D 2hc b#! = (h0c − k0c) exp " a D 2hc b# , (62)
where a ≥ 0 and b ≥ 0 are two material constants, and
D = −(s a ij + scij) ˙eij q 2 3˙emn˙emn + v u u u t (sa ij + scij) ˙eij q 2 3˙emn˙emn 2 + (ha+ hc)2− 3 2(s a ij + scij)(saij + scij) (63)
is the distance between stress point sij and its image point on bounding surface along the strain
rate direction ˙eij. The bounding surface is a surface in stress space with center sdij and radius
ha+ hc, whereas the contact surface is a surface with center sdij and radius hc. When fa = 0 and
fc = 0, the yield surface is in contact with the bounding surface.
The projection of saij + scij onto the unit tensor along the direction ˙eij is denoted by
d1 = (sa ij + scij) ˙eij q 2 3˙emn˙emn . (64)
Thus we obtain the following relation:
(D + d1)2+ d22 = (ha+ hc)2, (65) where d22 = 3 2(s a ij + s c ij)(s a ij + s c ij) − d 2 1. (66)
From Eqs. (64)-(66) we can derive Eq. (63). In addition, we suppose that
G > 0, ha> 0, hc > 0, h0c > 0, k 0
In this formulation a contact surface is employed to confine the motion of contact stress. While the on-off switching criteria of plasticity are derived from the first complementary trio, the switching criteria of kinematic hardening rules are derived according to the second complementary trio. In terms of the new concept of contact stress and contact surface, it becomes easy to derive the governing rule of back stress during the contact of yield surface and bounding surface. The validity and accuracy of the new model are confirmed by comparing the numerical results with the experimental data for SAE 4340 and RHA materials under four uniaxial cyclic testings and three biaxial cyclic testings. Even the material constants used in the model are parsimonious with only twelve, it is immediately recognized that the cyclic response curves described by the new model are in good agreement with experimental data. A detailed description may be found in [5].
8 Distributed switching mechanism
Mimicking the complementary trio continuum (33), we may devise the distributed switching mechanism as follows:
f (ν) ≤ 0, ˙λ(ν) ≥ 0, f(ν) ˙λ(ν) = 0 ∀ν ∈ I ⊆ R. (68) The distributed switching mechanism is versatile in modeling complex phenomena in plasticity and friction vibration problems [8].
9 Yield thick-surface
The multiple switching mechanism (34) and the distributd switching mechanism (68) may be employed to overcome the “over-square” problem at the elastic-plastic transition raised in section 1. However, the N -fold switching mechanism needs N yield surfaces fν = 0, ν = 1, · · · , N and the
distributed switching mechanism needs uncountably infinitely many yield surfaces f (ν) = 0, ν ∈ I ⊆ R. Therefore, as an alternative, we propose a complementary trio with a yield thick-surface:
hf i ˙f ≤ (1 − f ) ˙λ, ˙λ ≥ 0, hfi ˙f ˙λ = (1 − f ) ˙λ2, (69) where without loss of generality f and λ are assumed to be dimensionless and hf i = (f + |f |)/2. Note that (69) preserves the structure of a complementary trio because (69) becomes
upon requiring
φ = hf i ˙f − (1 − f ) ˙λ. (71)
It is assumed that {σ|f (σ, p, λ, β) ≤ 0} is simply connected in stress space, so is {σ|f (σ, p, λ, β) ≤
1}, and that {σ|f (σ, p, λ, β) ≤ 0} ⊂ {σ|f (σ, p, λ, β) ≤ 1}. Here σ denotes stress.
If there exists at t → −∞ a zero-value state when all (external and internal) state variables vanish, which is satisfied for essentially all cases, the complementary trio (69) automatically im-plies several novel features. First, f ≤ 1 is the admissible domain and f > 1 is the forbidden; more precisely {σ|f ≤ 1} is the stress admissible domain and {σ|f > 1} is the stress forbidden in stress space. Second, there is a yield thick-surface in stress space bounded by the outer surface {σ|f = 1} and the inner surface {σ|f = 0}. The yield thick-surface {σ|0 ≤ f ≤ 1} is a closed re-gion of nonzero volume (nonzero measure) in stress space. Recall that a conventional yield surface is of zero volume in stress space. Third, the consistency condition becomes f ˙f = (1 − f ) ˙λ.
10 A model of elastoplasticity with yield thick-surface
Let us consider a (comprehensive) constitutive model of local, rate-independent, small defor-mation, flow elastoplasticity with yield thick-surface as postulated below. The strain and stress σ are related by = e+ p, (72) e = ∂W c(σ, p, λ, β) ∂σ , (73) ˙p = a(σ, p, λ, β) ˙λ, (74) ˙ β = b(σ, p, λ, β) ˙λ, (75) hf i ˙f ≤ (1 − f ) ˙λ, (76) ˙λ ≥ 0, (77) hf i ˙f ˙λ = (1 − f ) ˙λ2. (78)
Here e, p, β and λ are called the elastic strain, plastic strain, internal variable and measure of
plasticity, respectively; in fact, they are all internal variables in a broad sense.
The model has a total of four prescribed material functions, namely the complementary stored energy Wc(σ, p, λ, β), plastic flow tensor field a(σ, p, λ, β), internal variable flow vector
and β, and have the properties: C(σ, p, λ, β) := ∂ 2Wc(σ, p, λ, β) ∂σ∂σ exists. (79) M(σ, p, λ, β) := [C(σ, p, λ, β)]−1 exists. (80) ∂f ∂σM ∂2Wc ∂σ∂λ+ ∂2Wc ∂σ∂βb + a > ∂f ∂λ + ∂f ∂βb . (81) {σ|f (σ, p, λ, β) ≤ 0} ⊂ {σ|f (σ, p, λ, β) ≤ 1}. (82)
Here C is called the elastic compliance and M the elastic modulus, both being symmetric fourth order tensors (with 21 independent components).
In the above , e, p and σ are symmetric second order tensors (with six independent
components); β is a vector (with nβ independent components); and λ is a scalar. All these second
order tensors, vector and scalar are functions of one and the same independent variable, which is taken either as the ordinary time or as the parameter (e.g. arc length) of the control path; however, for convenience, the independent variable, no matter what it is, will be simply called “time” and given the symbol t. A superimposed dot denotes differentiation with respect to time, d/dt.
In this section we use no signs or symbols for tensor products and scalar products, for there are no essential difficulties to distinguish them in context. For example,
e = ∂W c(σ, p, λ, β) ∂σ means e ij = ∂Wc(σ, p, λ, β) ∂σij ; ˙e = C ˙σ + · · · means ˙eij = Cijkl˙σkl+ · · · ; ∂f ∂σM ∂2Wc ∂σ∂βb = ∂f ∂σijM ijkl ∂2Wc ∂σkl∂βmb m.
By a symmetric second order tensor σ we mean σij = σji, for example, and by a symmetric fourth
order tensor M we mean Mijkl = Mklij = Mjikl= Mijlk, for example.
Without loss of generality, it is also tacitly postulated that with the above differential model, there is a time instant designated as t0, called the zero-value (or annealed or remoulded) time,
before and at which the material is in the zero-value (or annealed or remoulded) state in the sense that the relevant values (t0), e(t0), p(t0), σ(t0), β(t0) and λ(t0) are all zero, and 0 ∈
{σ|f (σ(t0), p(t0), λ(t0), β(t0)) ≤ 0}. In general the zero-value time t0, the initial time ti and the
11 Transcritical bifurcation
Transcritical bifurcation is a bifurcation with one codimension and is charateristic of exchange of stability. Its (expanded) normal form is
˙x = α1ux − α2x2 (83)
where x is state and u is control. Usually the normal form is referred to (83) with α1 = α2 = 1.
We extend them to α1 > 0, α2 > 0 (and even to α1 ≥ 0, α2 ≥ 0). Examining the x-u bifurcation
diagram of (83), we find that its fixed points are described by 0 = α1ux − α2x2 and stable fixed
points by the complementary trio:
α1u ≤ α2x, x ≥ 0, α1ux = α2x2. (84)
Let us identify x, u, α1, and α2 with ˙λ, ˙f , hf i, and (1 − f ), respectively. We then has a
non-equilibrium switching mechanism: ¨
λ = hf i ˙f ˙λ − (1 − f ) ˙λ2. (85)
Its corresponding equilibrium states are
hf i ˙f ≤ (1 − f ) ˙λ, ˙λ ≥ 0, hfi ˙f ˙λ = (1 − f ) ˙λ2. (86) Note that (86) is nothing but (69). It is rate independent. On the contrary, the rate dependent equation (85) acts like a smeared switching mechanism containing a viscous effect.
12 A different category of models
For the N -fold switching mechanism (34), if fν = 0 for a particular ν ∈ Ia ⊆ I = {1, · · · , N },
it is said that the ν-th yield condition is satisfied, critical, or potentially active, and it is readily deduced that
˙
fν ≤ 0 if fν = 0, (87)
˙
fν˙λν = 0 if fν = 0, (88)
or, written in vector form,
˙fa≤ 0, (89)
Proof: From (34)1, fν(t + dt) ≤ 0, that is, fν(t) + ˙fν(t)dt ≤ 0. Hence, if fν(t) = 0, then
˙
fν(t) ≤ 0; thus (87) is proved. Next, from (34)3, fν(t + dt) ˙λ(t + dt) = 0, that is, [fν(t) +
˙
fν(t)dt][ ˙λν(t) + ¨λν(t)dt] = 0. Hence, if fν(t) = 0, then ˙fν(t) ˙λν(t) = 0; thus (88) is proved.
It cannot be overemphasized that (87)-(90) are merely a consequence of (34) and that it would be redundant to put together as in some literature (e.g. [1]) the consequent along with the original,
f ≤ 0, λ ≥ 0,˙ f · ˙λ = 0, ˙fa≤ 0, ˙fa· ˙λa = 0,
in formulating a model with switching mechanism. But it is rather easy to see that it would turn out to be wrong if one (e.g. [6]) wrote
f ≤ 0, λ ≥ 0,˙ f · ˙λ = 0, ˙f ≤ 0, ˙f · ˙λ = 0,
because the stress level fν was limited inappropriately by ˙fν ≤ 0 when fν < 0.
However, if we requre
˙f ≤ 0, λ ≥ 0,˙ ˙f · ˙λ = 0, (91)
instead of
f ≤ 0, λ ≥ 0,˙ f · ˙λ = 0, (92)
then we have a different category of models for which {σ|f ≤ 0} is the stress admissible domain that is everywhere potentially active (on) and also potentially passive (off). In other words, the stress admissible domain is everywhere nothing but the yield thick-surface. Indeed, in such category fν = 0 is a limit surface rather than a yield surface. Note that (91) preserves the structure of
complementary trios, as (91) becomes
φ ≤ 0, λ ≥ 0,˙ φ · ˙λ = 0, (93)
upon requiring
φ = ˙f. (94)
A notable model belonging to this category is given by Sewell [9] as follows:
˙q = C ˙Q + N ˙λ, (95)
˙f = NTQ − H ˙˙ λ ≤ 0, λ ≥ 0,˙ ˙f ˙λ = 0. (96)
It is observed that this model can be obtained by replacing f in the model (43)-(44) by ˙f. As just indicated, the model (95)-(96) and the model (43)-(44) belong to different categories.
13 Equations of state and evolution of elastoplasticity
We now turn our attention from the yield surface, yield thick-surface, and switching mecha-nism to more broad areas including the elastic equation and plastic flow and hardening rules.
To begin with, a general form of elastoplastic model may be written as follows: q = ∂W c(Q, qp, λ, β) ∂Q + q p, (97) ˙qp = a(Q, qp, λ, β) ˙λ, (98) ˙ β = b(Q, qp, λ, β) ˙λ, (99) φ ≤ 0, λ ≥ 0,˙ φ · ˙λ = 0. (100) Here Q = (Q1, Q2, · · · , Qn), q = (q1, q2, · · · , qn), and qp = (q p 1, q p
2, · · · , qnp) denote the n-vectors
of generalized stress, generalized strain, and generalized plastic strain, respectively. The vector β consists of nβ internal variables. In general nβ 6= n. The scalar Wc is the complementary
stored energy. An overdot marks rate, namely derivative with respect to ordering, not necessarily physical, time t. All Q, q, qp, β, φ
ν, fν, λν are functions of time. The matrix-valued functions a
and b specify the plastic flow and the internal variable flow, respectively. The vectors φ, f and λ each have N components. The complementary trios (100) furnish the model with 2N alternative
phases, and for φ there are three categories as elaborated above for selection:
φν = fν(Q, qp, λ, β), (101)
φν = hfνi ˙fν− (1 − fν) ˙λν, (102)
φν = ˙fν. (103)
All components may belong to one category or components may be divided into different categories. In addition to the above three categories for rate-independent models, we may choose to replace the complementary trios (100) with
¨
λν = hfνi ˙fν˙λν− (1 − fν) ˙λ2ν. (104)
resulting in non-equilibrium smeared switching mechanism for rate-dependent models.
The above general form (97)-(100) has the flavor familiar in the literature of plasticity, while in the literature of dynamical systems there is a general form of expressions which looks somewhat similar:
y = Y(t, x, u), (106) where x, u, y are the vectors of state, input, and ouput, respectively. Eqs.(105)-(106) are usually called a state-space representation of a dynamical system. Referring to the constitutive equations (97)-(100), we may replace t in the system equations (105)-(106) by λ and append to the system equations a suitable set of complementary trios, resulting in the following state-space representation of a dynamical system with 2N alternative phases:
dx = V(λ, x, u)dλ, (107)
y = Y(λ, x, u), (108)
φ ≤ 0, dλ ≥ 0, φ · dλ = 0. (109)
Comparison shows that if the generalized stress Q is taken as input u, the generalized strain q as output y, and qp and β together as state x, then the general constitutive model (97)-(100) is merely a special case of the dynamical system (107)-(109). The modified state equation (107) describes the combined effects of plastic flow and all kinds of hardening and softening, while the modified output equation (108) describes not only the elastic effect but also the composition (and even coupling) of the elastic and plastic effects. These two equations serve as equations of evolution and state in a very general way.
14 Concluding remarks
We studied the theory of mechenisms of switching for describing the phenomena of switching between plastic and elastic phases observed in plasticity experiments and switching between sliding and sticking phases observed in friction experiments. The first focus was on the complementary trio and the corresponding switching mechanism. We reviewed several feasible approaches to the related topics.
In order to deal with the problems that the models with a conventional yield surface tends to predict a loop of stress-strain curve which is over square in reverse loading and whose back stress is over oscillatory in cyclic loading, we generalized the switching mechanism to a multiple swithching mechanism and to a distributed switching mechanism.
As an alternative to the multiple and distributed switching mechanisms, we proposed to thicken the yield surface so that the plastic-phase constitutive equations can develop their functions in more ample space in stress space. To remedy the loss of the consistency condition due to the
thickening, the yield thick-surface must be associated with a switching mechanism different from that associated with the conventional zero-thickness yield surface.
In this report we also generalized the state and output equations in the theory of dynamical control systems to the elastoplastic equations of evolution and state in a very general way.
References
[1] A. Borkowski, Analysis of Skeletal Structural Systems in the Elastic and Elastic-Plastic Range, Elsevier/PWN, Amsterdam/Warsaw, 1988.
[2] A. Fischer, An NCP-function and its use for the solution of complementarity problems, in D.-Z. Du, L. Qi, and R. S. Womersley (eds.), Recent Advances in Nonsmooth Optimization, World Scientific, Singapore, pp.88-105, 1995.
[3] H.-K. Hong and C.-S. Liu, Coulomb friction oscillator: modelling and responses to harmonic loads and base excitations, Journal of Sound and Vibration, Vol.229, pp.171-1192, 2000.
[4] H.-K. Hong and C.-S. Liu, Internal symmetry in the constitutive model of perfect elastoplasticity, International Journal of Non-Linear Mechanics, Vol.35, pp.447-466, 2000.
[5] H.-K. Hong, C.-S. Liu, and Y.-P. Shiao, Two complementary trios model and experimental simu-lations of SAE 4340 and RHA, submitted to International Journal of Plasticity.
[6] G. Maier, A matrix structural theory of piecewise-linear plasticity with interacting yield planes, Meccanica, Vol.5(1), pp.55-66, 1970.
[7] O. L. Mangasarian, Equivalence of the complementarity problem to a system of nonlinear equations, SIAM Journal on Applied Mathematics, Vol.11, pp.89-92, 1976.
[8] V. Palmov, Vibrations of Elasto-Plastic Bodies, Springer, Berlin, 1998.
[9] M. J. Sewell, Maximum and Minimum Principles, Cambridge University Press, Cambridge, Eng-land, 1987.
In the course of the research we made several attempts at formulating the so-called yield thick-surface. We record in the following two attempts.
1. We modify the complementary trio (8) to be
a ≥ 0, b ≥ 0, b ≤ 1, (110)
and correspondingly propose a so-called yield thick-surface and related switching mechanism in the following:
˙λ ≥ 0, Φ ≤ 0, ˙λφ ≥ 0, ˙λ ˙f Φ = 0, f H(φ) ≤ 0,˙ (111)
with
{σ|φ(σ, p, λ, β) ≤ 0} ⊆ {σ|Φ(σ, p, λ, β) ≤ 0}.
Here σ denotes stress and H(φ) is Heaviside’s unit step function.
Note that introduced here is a yield thick-surface in stress space bounded by the maximum yield surface {σ|Φ = 0} and the minimum yield surface {σ|φ = 0}. The yield thick-surface {σ|φ ≥ 0 ≥ Φ} is a closed region of nonzero volume (nonzero measure) in stress space. By setting φ = f = Φ the yield thick-surface reduces to a yield surface in the conventional sense, which is of zero volume in stress space, and hence the proposed (111) recovers to a conventional one (28).
Equations (111) indicate that Φ ≤ 0 is the admissible domain and Φ > 0 is the forbidden, and more precisely {σ|Φ ≤ 0} is the stress admissible domain and {σ|Φ > 0} is the stress forbidden in stress space.
2. Without loss of generality let f and λ be dimensionless. ˙λ ≥ 0, f ˙λ ≥ 0, f H(f ) ≤ (1 − f ) ˙λ,˙ f ˙λ = (1 − f ) ˙λ˙ 2.
Automatically implied here is that f ≤ 1 if there exists at t → −∞ a zero-value state when all (external and internal) state variables vanish.
The two attempts have a common drawback. They lose the structure of the complementary trio (8). Notice that (110) is no longer a complementary trio.
Incremental finite element formulation preserving elastoplastic switching
formalism
Hong-Ki Hong and Shin-Jang Sung Department of Civil Engineering
National Taiwan University Taipei, Taiwan Abstract
This paper studies elastoplastic analysis of a two- or three-dimensional body. The focus is on deriving an incremental formulation and finite element algorithm capable of accomodating con-stitutive models with kinematic hardening-softening, multiple yield modes, yield thick-surface, etc. The finite element algorithm obtained preserves the formalism of complementary trios that the original constitutive models have. An illustrative example is given to demonstrate the usability of the algorithm.
1 Introduction
The numerical solutions of plasticity problems cost a lot of computational efforts, since the nature of plastic deformation is progressive, alternative, and nonlinear. Due to the progressive nature, numerical calculations are carried out step by step at each point of the body; due to the alternative nature, monitoring of elastic or plastic phases is necessary at each step and each point; and owing to nonlinearity, in order to reach an acceptable accuracy, an iterative procedure is often invoked at each step and each point. To avoid iterations, Zhong et. al [3, 4] proposed an algorithm on the basis of the parametric variational principle. Their finite element formulations for plasticity problems preserve the so-called complementarity formalism. This feature was not at all typical in the literature of plasticity but did appear earlier in the papers of Maier et. al [1, 2] on “piecewise linear” plasticity. In the present paper we want to extend the formulation to allow for sophisticated plastic ingredients such as kinematic hardening-softening, multiple yield modes, yield thick-surface, etc. in a more versatile class of constitutive models. Moreover, instead of the parametric variational principle used in [3, 4], a simpler yet more logical derivation is carried out, especially for the complementary trios.
The class of constitutive models considered here is dσij = Eijkl(dε kl− dεpij), dεpij = aνij(σ, εp, λ, β) dλν, dβI = bνI(σ, εp, λ, β) dλν, fν(σ, εp, λ, β) ≤ 0, dλ ν ≥ 0, fνdλν = 0. (1) Here σij, ε
ij, εpij, Eijkl, fν are stress, strain, plastic strain, elastic modulus and yield functions,
respectively. Generally, i, j, k, l = 1, 2, 3 for problems of 3-dimensional domain D but i, j, k, l = 1, 2 for 2-dimensional problems, and ν, η = 1, 2, ..., Ny where Ny is the number of yield surfaces.
The I of βI (and bνI) is an index set of subscripts and superscripts, indicating that βI may be the
first order, the second order, the third order tensors, or the set of them. When these subscripts and superscripts appear in pairs in the same terms, the summation convention follows. In the models all εpij, λν, βI are internal variables. Material properties are specified by Eijkl, aνij, bνI, fν,
which may in general be taken as functions of the internal variables εpij, λν, βI as well as stress
σij; however, in many applications, Eijkl reduces to be a constant tensor, and in some other
specific applications aν
ij can be the gradient ∂gν/∂σij of the plastic flow potential functions
gν(σ, εp, λ, β), or even the gradient ∂fν/∂σij of the yield functions fν.
Note that the models (1) include kinematic and isotropic hardening rules, multiple yield surfaces, and, moreover, yield thick-surfaces if we
replace fν by hfνidfν − (1 − fν)dλν (no sum on ν) (2)
where hfνi = (fν + |fν|)/2. Without introducing (2) the yield surfaces are fν = 0, of zero
thicknesses. On the contrary, it is with (2) such that the yield surfaces are the regions between fν = 0 and fν = 1, the so-called yield thick-surfaces.
For the elastoplastic models (1) the stresses corresponding to the elastic strain εe
ij = εij−εpij
and to the plastic strain εpij are the same, as can be seen from (1)1. Although the
mechanical-element model shows that there are two mechanisms apart, a material point in a body has two mechanisms simultaneously. That implies the stress distributions in both of them are the same, too. However, their responses are additive. Just like a series-of-spring system, the principle of equilibrium indicates that every spring indeed is subjected to the same force and that, due to the prescribed relation of kinematics, their responses are additive. Therefore, in the following formulation we divide our derivation into two parts, the part corresponding to the elastic strain and the part corresponding to the plastic strain. In the part correponding to the elastic strain we convert the strong form of equilibrium equations with boundary conditions into an incremental weak form suitable for finite element discretization and incorporate the elastic
relations. As for the part corresponding to the plastic strain, we reform the yield functions and the complementary trios into an incremental weak form and incorporate the plastic relations.
On the basis of this point of view, the incremental weak formulation is derived in Section 2. With some assumptions a finite element algorithm is presented in Section 3. An illustrative example is worked out in Section 4.
2 Incremental variational formulation
The variables of the elastoplastic models are path dependent. Let us denote the path pa-rameter by t and call it time for simplicity, although it is an ordering scalar for the purposes of constitutive models of elastoplasticity, not necessarily (an ordinary physical) time, but it can always be chosen to be (an ordinary physical) time. As such σij, ε
ij, εpij, βI, fν, λν, etc. at time
t are respectively expressed as σij(t), ε
ij(t), εpij(t), βI(t), fν(t), λν(t), etc. Now we study a body
with material (1) subjected to the body force bi prescribed in the domain D and to the traction ¯
pi and displacement ¯ui prescribed on the boundary B. The prescribed bi, ¯pi, and ¯ui are not only
functions of spatial position x but also functions of (the ordinary physical) time t. Therefore, it is obvious all the variables mentioned above and to be named below must be functions of (the ordinary physical) time t and spatial position x. However, the dependence on x will not be explicitly written out in the context, except for (22), for the sake of brevity. Let ui(t) and
¨
ui(t) denote the displacement and acceleration, repectively, at time t. Step by step calculations
demand three kinds of quantities: at time t, at time t + ∆t, and in ∆t. Defining incremental variables as the differences of the variables evaluated at between time t and t + ∆t, we have ∆σij(t), ∆εij(t), ∆ε
p
ij(t), ∆βI(t), ∆fν(t), ∆λν(t), ∆bi(t), ∆¯pi(t), ∆¯ui(t), ∆ui(t), ∆¨ui(t), etc.
Let us first consider the part corresponding to the elastic strain. For the kinetic variables σij(t) and the kinematic variables ui(t), the equations of equilibrium in D and the boundary
conditions on Bp and Bu at time t are
−σij,j(t) − bi(t) + ρ¨ui(t) = 0 in D, σij(t)n j− ¯pi(t) = 0 on Bp, ui(t) − ¯ui(t) = 0 on Bu, (3)
respectively, where ρ is the mass density of the body. To prepare for the finite element algorithm to be given in the next section, modifications of (3) are needed. At first, the incremental version of (3) is −∆σ,jij(t) − ∆bi(t) + ρ∆¨ui(t) = 0 in D, ∆σij(t)nj− ∆¯pi(t) = 0 on Bp, ∆ui(t) − ∆¯ui(t) = 0 on Bu. (4) 3
By choosing an appropriate ∆ui(t), (4)3 can be automatically satisfied. With the aid of
essen-tially arbitrary functions ri(t), often called testing functions, the strong form (4), being for a
typical material point, may be rewritten equivalently in integral (weak) form: Z D ri(t)[−∆σ,jij(t) − ∆b i(t) + ρ∆¨ui(t)]dD + Z Bp ri(t)[∆σij(t)nj − ∆¯pi(t)]dB = 0. (5)
Then, by letting ri(t) = δ∆ui(t) and according to the divergence theorem
Z D ∆σ,jij(t)δ∆ui(t)dD = − Z D ∆σij(t)δ∆ui,j(t)dD + Z B (∆σij(t)nj)δ∆ui(t)dB, (6)
the integral equation (5) becomes Z
D
[δ∆ui,j(t)∆σij(t) + δ∆ui(t)(−∆bi(t) + ρ∆¨ui(t))]dD
− Z Bu δ∆ui(t)∆σij(t)njdB − Z Bp δ∆ui(t)∆¯pi(t)dB = 0. (7)
However, δ∆ui(t) must vanish on Bu according to (3)3. Hence, the above equation (7) combined
with (1)1 and (1)2 can be expressed as
Z
D
[δ∆εij(t)Eijkl(∆εkl(t) − aklν∆λν(t)) + δ∆ui(t)ρ∆¨ui(t)]dD
= Z Bp δ∆ui(t)∆¯pi(t)dB + Z D δ∆ui(t)∆bi(t)dD. (8)
It will be seen later that ∆λν(t) serves to be a connection relating the elastic equations (8) to
the plastic expressions to be elaborated right below.
Next, let us turn our focus to the part corresponding to the plastic strain. With the notations of the incremental variables, the complmentary trios (1)4 become
fν(t + ∆t) ≤ 0, ∆λν(t) ≥ 0, fν(t + ∆t)∆λν(t) = 0. (9)
for the time period between t and t + ∆t, where
fν(t + ∆t) = fν(t) + ∆fν(t) = fν(t) + Wνkl(t)∆εkl(t) + Zνη(t)∆λη(t) (10) in which Wνkl := ∂f ν ∂σijE ijkl, (11) Zνη := −∂f ν ∂σijE ijklaη kl+ ∂fν ∂εpija η ij + ∂fν ∂λη +∂f ν ∂βI bηI. (12) 4
In the above derivation use has been made of the elastic and plastic relations (1)1-(1)3.
Multiplying fν and ∆λν respectively by positive constants cf and cλ with reciprocal
di-mensions such that cffν and cλ∆λν are nondimensionalized and so are (9), giving
cffν(t + ∆t) ≤ 0, cλ∆λν(t) ≥ 0, cffν(t + ∆t)cλ∆λν(t) = 0. (13)
It can be readily shown that the dimensionless complementary trios (13) are equivalent to cλ∆λν(t) − cffν(t + ∆t) − |cλ∆λν(t) + cffν(t + ∆t)| = 0. (14)
With the aid of nonnegative testing functions rν(t), which are nonnegative but otherwise
arbi-trary and mutually independent, the strong form (14) is equivalent to the weak form: Z
D
[ rν(t)cλ∆λν(t) − rν(t)cffν(t + ∆t) − |rν(t)cλ∆λν(t) + rν(t)cffν(t + ∆t)| ] dD = 0. (15)
From (13) and the nonnegativity of rν(t) we may divide the domain D into De(t) and Dp(t)
such that rν(t)cλ∆λν(t) > 0 and rν(t)cffν(t + ∆t) = 0 in Dp(t), (16) rν(t)cλ∆λν(t) = 0 and rν(t)cffν(t + ∆t) ≤ 0 in De(t). (17) In Dp (16) ensures that Z Dp [ rν(t)cλ∆λν(t) − rν(t)cffν(t + ∆t) − |rν(t)cλ∆λν(t) + rν(t)cffν(t + ∆t)| ] dD = Z Dp rν(t)cλ∆λν(t)dD − | Z Dp rν(t)cλ∆λν(t)dD|, (18) Z Dp rν(t)cffν(t + ∆t)dD = 0. In De (17) ensures that Z De [ rν(t)cλ∆λν(t) − rν(t)cffν(t + ∆t) − |rν(t)cλ∆λν(t) + rν(t)cffν(t + ∆t)| ] dD = − Z De rν(t)cffν(t + ∆t)dD − | Z De rν(t)cffν(t + ∆t)dD|, (19) Z De rν(t)cλ∆λν(t)dD = 0. 5
Substituting (18) and (19) into (15) we have Z D rν(t)cλ∆λν(t)dD− Z D rν(t)cffν(t+∆t)dD−| Z D rν(t)cλ∆λν(t)dD+ Z D rν(t)cffν(t+∆t)dD| = 0.
This is equivalent to the complementary trios: Z D rν(t)cffν(t + ∆t)dD ≤ 0, Z D rν(t)cλ∆λν(t)dD ≥ 0, Z D rν(t)cffν(t + ∆t)dD Z D rν(t)cλ∆λν(t)dD = 0.
Deletion of the positive constants gives Z D rν(t)fν(t + ∆t)dD ≤ 0, Z D rν(t)∆λν(t)dD ≥ 0, Z D rν(t)fν(t + ∆t)dD Z D rν(t)∆λν(t)dD = 0.
Letting rν(t) = δ∆λν(t) ≥ 0 we have the variational form
Z D δ∆λν(t)fν(t + ∆t)dD ≤ 0, Z D δ∆λν(t)∆λν(t)dD ≥ 0, Z D δ∆λν(t)fν(t + ∆t)dD Z D δ∆λν(t)∆λν(t)dD = 0. (20)
Substitution of (10) into the left-hand side of (20)1 gives
Z D δ∆λν(t)fν(t + ∆t)dD = Z D δ∆λν(t)[fν(t) + Wνkl(t)∆εkl(t) + Zνη(t)∆λη(t)]dD. (21)
Now we have available an incremental variational form, namely (8) and (20) with (21), suitable for finite element discretization.
3 Finite element algorithm
To take advantage of the finite element method (e.g. [5]) let
∆u = ∆u(t, x) = N(x)∆U(t) = N∆U, ∆ε = ∆ε(t, x) = B(x)∆U(t) = B∆U, ∆λ = ∆λ(t, x) = L(x)A∆Λ(t) = LA∆Λ, f0 = f (t, x) = L(x)AF(t) = LAF0, f = f (t + ∆t, x) = L(x)AF(t + ∆t) = LAF. (22) 6
Here the capital symbols ∆U, ∆Λ, F0, and F denote column matrices of nodal values, N a
matrix of interpolation functions for ∆ui and shape functions of finite elements, B the
strain-displacement matrix, L a matrix of interpolation functions for ∆λν and fν. Inserted in (22)
after each L(x) is a constant square matrix A, which will be determined later on. With (22) we first discretize (8) as follows:
Z
D
[Nδ∆U]Tρ∆ ¨UdD + Z
D
[Bδ∆U]TE(−aLA∆Λ + B∆U)dD
= Z Bp [Nδ∆U]T∆¯pdB + Z D [Nδ∆U]T∆bdD.
Because the virtual incremental nodal displacement δ∆U is arbitrary and independent of x, we can drop out δ∆U so that the above equation becomes
M∆ ¨U − V∆Λ + K∆U = ∆P, (23) upon defining M = X Z De NTρdD, V = X Z De BTEaLAdD, K = X Z De BTEBdD, ∆P = X Z Bpe NT∆¯pdB +X Z De NT∆bdD.
The subscript “e” means that these symbols describe only one element, not to confuse it with the elastic part. The symbol P does not mean summation but assembly; thus, for a finite element occupying a domain De during a time period between t and t + ∆t, first integrate over
each and every finite element and then assemble. Next, we discretize (20) as follows:
Z D [LAδ∆Λ]TLAFdD ≤ 0, Z D [LAδ∆Λ]TLA∆ΛdD ≥ 0, Z D [LAδ∆Λ]T LAF TdD Z D [LAδ∆Λ]TLA∆ΛdD = 0. (24) 7
Similarly, because δ∆Λ is arbitrary and independent of x, let us drop out δ∆Λ so that Z D [LA]TLAFdD ≤ 0, Z D [LA]TLA∆ΛdD ≥ 0, Z D [LA]TLAF T dD Z D [LA]TLA∆ΛdD = 0. (25)
This can be written as
F ≤ 0, ∆Λ ≥ 0, FT∆Λ = 0, (26) upon setting A = AT = XZ De LTLdD −1/2 . (27)
To discretize (21) we follow the same procedure as from (24) to (27), obtaining
F = F0+ G∆U + H∆Λ, (28) where G = X Z De [LA]TLA(∂f ∂σ) TEBdD, H = X Z De [LA]TLA −(∂f ∂σ) TEa + (∂f ∂εp)a + ∂f ∂λ + ( ∂f ∂β) Tb LAdD.
Now, collecting (23), (26), and (28), we obtain the finite element formulation: M∆ ¨U − V∆Λ + K∆U = ∆P, F = F0+ G∆U + H∆Λ, F ≤ 0, ∆Λ ≥ 0, FT∆Λ = 0. (29)
Comparison shows that the finite element formulation (29) preserves the formalism of the complementary trios the constitutive models (1) have.
Because the response displacement and strain can be divided into the elastic and plastic parts which are additive, we may first set ∆Λ ≥ 0 to get the elastic displacement from the first equation in (29). However, these responses, a middle product of finding stress distribution, is not real since the plastic responses are not yet considered. After finding stress distribution, G and H can be determined so that there are only two column matrices ∆U and ∆Λ left unknown. Then, by substituting the first equation of (29) into the second one, the unknown ∆U is cancelled and
F = F0+ GK−1∆P + [GK−1V + H]∆Λ. (30)
Again assume ∆Λ = 0 in (30) to calculate F. If some entries in F can not satisfy the inequalities in (29), it means some elements are in the plastic phase. In that case, set the illegal entries of F to 0 and use the whole (30) to determine ∆Λ. However, it is of concern that in order to reduce degrees of freedom, the rows in (30) with values of F legal should not participate in the calculation. The response of the time step is then determined by summing up the elastic and plastic parts. The dimensions of the rectangular matrices are summarized in Table 1.
number of number of row column N Nd Nf ∗ Ep B Nd(Nd+1) 2 Nf ∗ Ep E Nd(Nd+1) 2 Nf(Nf+1) 2 L Ny Ey A Ey Ey Me Nf ∗ Ep Nf ∗ Ep Ve Nf ∗ Ep Ny ∗ Ey Ke Nf ∗ Ep Nf ∗ Ep Ge Ny ∗ Ey Nf ∗ Ep He Ny ∗ Ey Ny ∗ Ey M Nf ∗ Nη Nf ∗ Nη V Nf ∗ Nη Ny ∗ Nλ K Nf ∗ Nη Nf ∗ Nη G Ny ∗ Nλ Nf ∗ Nη H Ny ∗ Nλ Ny ∗ Nλ symbol definition
Eη number of elements of the problem
Nη number of u-nodes of the problem
Ep number of u-nodes of each element
Nd degrees of freedom of displacement field
Nf degrees of freedom of each u-node
Nλ number of f-nodes of the problem
Ey number of f-nodes of each element
Ny number of yield surfaces
Table 1: matrix dimension
4 Example
In this section the finite element algorithm obtained is applied to an example of a 2-dimensional cantilever plate for illustration purposes. To simplify, we assume the problem to be static and the material has only one yield surface - von Mises yield condition; thus ν = 1 only, so the index ν can be dropped out. We choose triangular elements for shape, linear
ments for ∆ui, and constant elements for ∆λ and f . Prager’s kinematic hardening rule is used
here and the stiffness after yielding is 1/10 of the elastic stiffness. The material properties are: Young’s modulus E = 70 GPa, Poisson ratio = 0.3 and the yield stress in pure shear k = 270 GPa. The geometry properties are: length L = 150 mm, height h = 50 mm and the left end is fixed. The distributed line load q = 40 MPa∗m is applied on the top of the plate.
The result is shown in Figure 1 where there are three types of shape: the original shape, the elastic displacement and the total displacement. Figure 2 shows which element is in the plastic phase and the darker elements have more plastic strains. The vertical elastic displacement of the right lowest point is −0.00341 m. A rough check is to use the Bernoulli beam theory and the result is −0.00347 m. The vertical total displacement of the right lowest point is −0.00688 m.
5 Concluding remarks
In this paper it was shown that a comprehensive class of elastoplastic constitutive models can be incorporated into an incremental variational formulation, based on which a finite element discretization and step by step calculations can be performed. The illustrative example demon-strated that the derivation and the finite element algorithm are workable. But the example has only one single yield surface of von Mises. When using a new constitutive model different from the model of the example, only simple modifications are needed to fit the demand if the new model belongs to the class of models (1) which is comprehensive in scope. This advantage is more obvious in computer programming. The finite element formulation obtained preserves the formalism of complementary trios and falls into the type of complementarity problems. The present paper did not delve into the topic of solution methods for the complementarity problems. In the example the solution method used was relatively primitive. It is suggested that the solution methods for the type of problems need more studies, especially for problems with many yield surfaces. Moreover, the methods are not the only important strategic point. We should also explore what applications the proposed formulation may be applied to and the potential of it will be discovered after that.
References
[1] Corradi, L., On compatible finite element models for elastic plastic analysis, Meccanica, Vol.13, pp.133-150, 1978.
[2] Maier, G. and Nappi, A., On the unified framework provided by mathematical program-ming to plasticity, in D. J. Dvorak and R. T. Shield (eds.) Mechanics of Material Behavior, Elsvier, Amsterdam, pp.253-273, 1984.
[3] Zhang, H. W., Xu, W. L., Di, S. L. and Thomson, P. F., Quadratic programming method in numerical simulation of metal forming process, Computer Methods in Applied Mechanics and Engineering, Vol.191, pp.5555-5578, 2002.
[4] Zhang, H. W., Zhang, X. W. and Chen, J. S., A new algorithm for numerical solution of dynamic elastic-plastic hardening and softening problems, Computers and Structures, Vol.81, pp.1739-1749, 2003.
[5] Zienkiewicz, O. C. and Taylor, R. L., The finite element method Volume 1 : The basis, Elsevier Butterworth-Heinemann, Burlington, 2000.
Figure 1: The original and displaced positions of finite element nodes. Also shown are the elastic displacements.
Figure 2: Dark elements are in the plastic phase and the darker elements have more plastic strains.
Estimation of Model Parameters Using Group Preserving
Scheme
Chein-Shan Liua Li-Wei Liub, Hong-Ki Hongb
a Department of Mechanical and Mechatronic Engineering, Taiwan Ocean University,
Keelung, Taiwan
b Department of Civil Engineering, Taiwan University, Taipei, Taiwan
Abstract
In this report, we propose several methods for estimation of model parameters. Using the concepts of evolutionary equation theory and eigenfunction expansion, we convert an infinite dimemsional system into a system of ordinary differential equations, which are then augmented into Lie algebra of Lorentz group. A group preserving scheme solves the ordinary differential equations obtaining a set of algebraic formulae, based on which parameters are estimated. The finite element method is also used to approx-imate the infinite dynamical system and a group preserving scheme is also developed to estimate parameters. For constitutive laws of elastoplasticity, we propose methods for identification on Minkowski spacetime, for which least squares and Kalman filter theory are adopted to estimate the parameters. For the boundary value problem of plasticity, an incremental variational formulation preserving elastoplastic switching for-malism which the constitutive model has are derived. On the basis of the formulations, the corresponding method of identification is developed. Numerical results are shown in good agreement with theoretical results.
1. Introduction
In order to solve the inverse heat conduction problems, several methods have been proposed in the last decade. Huang and his coauthors dealt with three-dimensional inverse heat conduction problems by using the conjugate gradient method [7, 8] and the steepest descent method [9, 10] in conjunction with commercial codes. Telejko and Malinowski [16] converted the partial differential equation of heat conduction into a system of algebraic equations by using the finite element method in spatial domain and Galerkin integration scheme in time domain. Then, the Broyden-Fletcher-Goldfarb-Shanno method is used to deal with the optimum problem of estimating the thermal conductivity of two-dimensional problems. Cabeza et. al. [1] used Duhamel integral to express the dependence between boundary flux and response of temperature and approximated it by the system of algebraic equations. They employed the sequential
singular value decomposition to overcome the ill-condition problem and determine the boundary flux. Girault and Petit [4, 5] converted the heat equation into a system of algebraic equations by using finite volume method in spatial domain and Euler implicit scheme in time domain. They proposed a model reduction method with iterative least square method to deal with the optimum problem of determining the boundary flux.
Without iteration, Del Barrio [3] borrowed the concept of optimal quadratic con-trol. They converted the partial differential equation into the system of differential equations by using finite element method and then solved the direct, adjoint and sen-sitivity problems of the optimum problem with Tikhonov regularization in the same time. Yeung and Lam [17] converted the partial differential equation into the algebraic equations directly by using finite difference method and estimated the unknown param-eters by solving these algebraic equations. Lin et.al. combine method of Yeung and Lam with least square concept to deal with the inverse heat conduction problems.
In the elastoplastic problem, Hong and Liu [11] derived the exact solutions of model of perfect elastoplasticity with a non-zero constant rate of generalized strain path. They also derived a group preserving scheme for general paths.
In this article, we propose a method of estimating the model parameters. In the second section, we convert the partial differential equation into the evolutionary equation. We introduce the group preserving scheme and used it to derive the algebraic relation between two states of different time.
In the third section, the finite dimensional approximation of evolutionary equation is proposed. Also, the algebraic relation between two states of different time is derived by using group preserving scheme.
On the finite dimensional approximation, the exact method of identification cor-responding to the algebraic relation is showed in the forth section. In the fifth section, we addressed the methods for identification on the Minkoeski spacetime.
For the boundary value problem of elastoplasticity, the method of identification is proposed according to the formulations which preserve the elastoplastic switching formalism in the sixth section. The numerical examples are showed in the final section.
The well-known heat equation is written as follows: ∂ ∂tu(x, t) = ∇ · (k(x)∇u(x, t)) + g(x, t), t > 0, x ∈ Ω, α(x)u(x, t) + β(x)∂x∂ u(x, t) = γ(x), t > 0, x ∈ ∂Ω, u(x, 0) = u0(x), x ∈ Ω. (1) If we define a operator A as Aϕ = ∇ · (k(x)∇ϕ), ϕ ∈ D(A) = {ϕ ∈ H2(Ω); αϕ(x) + βϕ0(x) = γ, x ∈ ∂Ω}, in which H2(Ω) denotes the space of functions in L2(Ω) whose distribution derivatives
of order ≤ 2 are in L2(Ω) and L2(Ω) denotes Hilbert space with inner product hϕ, ψi := (RΩϕ(x)ψ(x)dx)1/2 and the norm kϕk := (hϕ, ϕi)1/2, then Eq.(1) can be written as follows:
ut= Au + g
u(0) = u0
(2) If the discrete eigenvalues λi of the operator A and the corresponding orthonormal
eigenvectors φi, i = 1, 2, · · · exist, then the u, g, and operator A can be written in term
of eigenvectors φi as follows: u =X i uiφi, g =X i giφi, Au =X i λiuiφi.
where ui := hu, φii and gi := hg, φii. Then the evolutionary equation can be converted
to the set of ordinary differential equations: d
dtui = λiui+ gi, i ∈ Z
+. (3)
According to Parseval equality,
hu, ui := Z Ω u2dx =X i u2i, we obtain d dtkuk = 1 (P iu2i)1/2 X (λiui+ gi)ui. (4)
Combining Eqs.(3) and (4), we can derive an augmented dynamical system ˙ X := d dt u kuk = O h hT 0 u kuk =: AX (5)
where u := u1 u2 .. . , kuk := ( X i u2i)1/2, h := λ1u1+g1 (P iu2i)1/2 λ2u2+g2 (P iu2i)1/2 .. . .
The state vector X always preserves the cone conditions
XTgX = 0, (6) where g = I 0 0 −1 ,
is a Minkowski metric, I is the identity matrix, and the state matrix A satisfies
ATg + gA = 0. (7)
The solution of the augmented dynamical system (5) can be expressed as follows: X(t) = G(t)G−1(t)X(t1), (8)
in which, G(t), known as the fundamental solution of (5) satisfying ˙
G(t) = A(t)G(t), (9)
G(0) = I. (10)
It is easy to show the fundamental solution G(t) preserves the Minkowski metric [11] GT(t)gG(t) = g. (11) Hence, the fundamental solution G(t) of Eq.(5) belongs to the Lorentz group. Besides, the A which satisfies condition (7) is a element of the Lie algebra [2]. For every element a of Lie algebra, there exist an associated element, exp(at), of Lie group satisfying
a = ( d
dt exp(at))|t=0. (12) The exp(at) is known as the exponential map of a. According to Eqs.(9), (10), and (12), the fundamental solution G(t) is a element of Lie group associated with the element A of Lie algebra.
Assuming that, the state matrix A is contant during t ∈ (t0, tT), the transition
map, G(tT)G−1(tT), can be derive
Hence, we obtain the algebraic equations between ui(t0) =: u0i and ui(tT) := u0i as follows: uTi = u0i + b( P i(u 0 i)2)1/2 (P i(λ 0 iu0i + g0i)2)1/2 (λ0iu0i + gi0) + (a − 1) P i(λ 0 iu0i + gi0)u0i P i(λ 0 iu0i + gi0)2 (λ0iu0i + g0i), (14) (X i (uTi )2)1/2 = b P i(λ0iu0i + g0i)u0i (P i(λ 0 iu0i + g0i)2)1/2 + a(X i (u0i)2)1/2, (15) where a := cosh((tT − t0)( P i(λ 0 iu0i + g0i)2)1/2 (P i(u0i)2)1/2 ), b := sinh((tT − t0)( P i(λ 0 iu0i + gi0)2)1/2 (P i(u0i)2)1/2 ). 3. Finite dimensional approximation of evolutionary equation
In the last section, we have converted the evolutionary equation to the infinite dynamical system. However, it is hard to find all the of eigenvectors of operator A practically. In this section, we will propose a finite dimensional approach to approximate the solution of equation of Eq.(1).
Suppose that DM(A) is a finite dimensional subspace of D(A). Let {ψi}Mi=1 be a
basis of DM(A). Multiplying the first equation of Eq.(1) by ψj and integrating it on Ω,
we can derive the weak form Z Ω ut(x)ψjdv = − Z Ω ∇ψj· (k(x)∇u(x))dV + Z Ω g(x, t)ψjdV + Z ∂Ω ψj(k(x)∇u(x)) · nds. (16) We project u onto DM(A) and write
u =
M
X
i=1
aiψi. (17)
Substituting Eq.(17) into Eq.(16) and dividing the domain Ω and boundary ∂Ω into M subdomains {Ωe}Me=1 and M portions of boundaries {∂Ωe}Me=1, respectively, we can
derive finite dimensional dynamical system as follows:
where H = (Hij) , Hij := − M X e=1 Z Ωe ∇ψi· (k(x)∇ψj)dV, ∀i, j = 1, · · · , M, (19) T = (Hij) , Tij := − M X e=1 Z Ωe ψiψjdV, ∀i, j = 1, · · · , M, (20) at := ha1 a2 · · · aMit (21) rt= hr1 · · · ri · · · rMit, ri := M X e=1 Z Ωe g(x, t)ψjdV − M X e=1 Z ∂Ωe ψi(k(x)∇u(x)) · nds, ∀i = 1, · · · , M. (22)
Consider the dynamical system, we can construct the state vector X and state matrix A of its M + 1-dimensional augemented dynamical system
d dtX(t) :=˙ d dt a kak = 0M ×M f kak ft kak 0 ! a kak =: AX(t). (23) satisfying the cone condition,
XT(t)gX(t) = 0, (24) and
ATg + gA = 0, (25)
where
f := T−1Ha + T−1r.
Hence, the fundamental solution G(t) of Eq.(23) belongs to the Lorentz group.
Similarly, the Eq.(23) possesses the same structure of Eq.(5). Assuming that, the state matrix A is constant during t ∈ (t0, tT), the transition map, G(tT)G−1(t0), can
be approximated as G(tT)G−1(t0) ∼= exp(A(tT − t0)) = I + a(tT,t0)−1 kf (t0)k2 f (t0)f T(t 0) b(tTkf (t,t00)f (t)k0) b(tT,t0)fT(t0) kf (t0)k a(tT, t0) ! , (26) where a(tT, t0) := cosh( (tT − t0)kf (tT, t0)k ka(t0)k ), b(tT, t0) := sinh( (tT − t0)kf (tT, t0)k ka(t0)k ). (27) Because the transition matrix satisfies condition (11), it preserves the Minkowski metric and belongs to the Lorentz group. Therefore, we obtain the algebraic equations between X(t0) and X(tT) X(tT) = I +a(tT,t0)−1 kf (t0)k2 f (t0)f T(t 0) b(tT,t0)f (t0) kf (t0)k b(tT,t0)fT(t0) kf (t0)k a(tT, t0) ! X(t0). (28)
4. Identification using group preserving scheme
In the last section, we have converted the dynamical system (18) to the algebraic equations (28). In this section, we will estimate the parameters of the dynamical system by using the algebraic equations.
Let η := a(tT, t0) − 1 kf (t0)k2 fT(t0)a(t0) + b(tT, t0)ka(t0)k kf (t0)k , (29)
we can rewrite the first n equations of Eq.(28) as follows:
a(tT) = a(t0) + ηf (t0). (30)
For identification, the states X(t0) and X(tT) are only known. According to
Eq.(29), we know that η is unknown because of the unknown f (t0). However, in the
group preserving scheme, there is another information on the last equation of Eq.(28) in addition to Eq.(30). Thereinafter, we will derive an exact formulation of η only in terms of x(t0) and x(tT).
According Eq.(30), we obtain f (t0) =
1
η(a(tT) − a(t0)).
Extracting the last equation of Eq.(28) and dividing it by kx(t0)k, we obtain
ka(tT)k ka(t0)k = b(tT, t0) fT(t 0)a(t0) kf (t0)kka(t0)k + a(tT, t0). (31) Let cos θ := f T(t 0)a(t0) kf (t0)kka(t0)k , s := (tT − t0)k(a(tT) − a(t0))k ka(t0)k , (32)
and from Eqs.(27) and (31) it follows that ka(tT)k ka(t0)k = cosh(s η) + cos θ sinh( s η). (33) Defining z := exp(s η),
we rewrite Eq.(33) to a quadratic equation of z as follows: (1 + cos θ)z2+ −2ka(tT)k
ka(t0)k
