Belytschko T. - Introduction (779635), страница 46
Текст из файла (страница 46)
(B4.5.11), the last term vanishes. The first term on the RHS can bereduced to the traction boundary since δui = 0 on Γu0i and Γt0i = Γ0 −Γu0i , sonSD∂0δuPdΩ=δunPdΓ=δui t i0dΓ0()∫Ω ∂X j i ji 0 Γ∫ i j ji 0 ∑∫i=1 Γ 000t(4.8.5)iwhere the last equality follows from the strong form. From Eq. (3.2.14)we note that ∂u δFij =δ i ∂X j (4.8.6)Substituting Eq. (4.8.5) into (4.8.3) and the result into (4.8.2) gives, after a change of sign andusing (4.8.6):n SD0∫ (δFij Pji −δui ρ 0bi + δui ρ0u˙˙ i )dΩ0 − ∑ ∫ δuiti dΓ0 = 0Ω0(4.8.7)i=1 Γ 0tiorn SD0∫ (δF :P − ρ0 δu⋅b + ρ0 δu⋅ ˙u˙ )dΩ0 − ∑ ∫ (δu⋅e i ) (e i ⋅ ti )dΓ0 = 0TΩ0(4.8.8)i=1 Γ 0tiThe above is the weak form of the momentum equation, traction boundary conditions, and interiorcontinuity conditions. It is called the principle of virtual work, since each of the terms in Eq.(4.8.7) is a virtual work increment.
The weak form is summarized in Box 4.6, in which physicalnames are ascribed to each of the terms.This weak form can also be developed by replacing the test velocity by a test displacement inEq. (4.3.9) and transforming each term to the reference configuration. The total Lagrangian weakform, Eq. (4.8.8), is thus simply a transformation of the updated Lagrangian weak form.4-63T. Belytschko, Lagrangian Meshes, December 16, 1998Box 4.6Weak Form for Total Lagrangian Formulation:Principle of Virtual WorkWEAK FORM: if u ∈U andδW int (δu, u) − δ W ext (δu,u ) +δ W inert (δu,u ) = 0∀δu∈U0(B4.6.1)then equilibrium, the traction boundary conditions and internal continuity conditionsare satisfied.In the above∫δW int = δFT : P dΩ0 =Ω0δW ext =∫ δFij Pji dΩ 0Ω0nSD(δu⋅e i )(e i ⋅ ti )dΓ0 = 0∫ ρ0δu ⋅bdΩ0 + ∑∫i =1Ω0=δW0Γt0inSD∫ δui ρ0 bi dΩ0 + ∑∫ δui t i dΓ0i =1(B4.6.3)= ∫ δu ⋅ρ0˙u˙ dΩ0 = ∫ δui ρ0˙u˙ idΩ0(B4.6.4)Ω0inert(B4.6.2)0Γt0iΩ0Ω0Weak Form to Strong Form.
Next we deduce the strong form from the weak form. To avoidwriting summations, we shall assume that on any part of the boundary, all traction or displacementcomponents are prescribed; the reader can easily generalize the proof to mixed boundary conditionswhere for each component, either the displacement or traction component is prescribed.Substituting (B4.8.6) into the first term of (B4.8.7) and the derivative product rule gives∫Ω0∂ (δui )P dΩ =∂Xj ji 0 ∂∂Pji δui Pji − δui dΩ∂Xj∂X j 0Ω0(∫)(4.8.9)Gauss's theorem on the first term on the RHS then yields:∫Ω0nSD∂Pji∂ (δui )Pji dΩ 0 = ∑ ∫ δui n0j Pji dΓ0 + ∫ δui n 0j Pji dΓ0 − ∫ δuidΩ0∂Xj∂X ji=1 Γ0Γ0Ωtiint0where the surface integral is changed to the traction boundary because δui = 0 on Γu0i andΓt0i = Γ 0 −Γu0i .Substituting (4.8.9) into (4.8.7) and collecting terms gives4-64T.
Belytschko, Lagrangian Meshes, December 16, 1998nSD ∂Pji˙˙0 = ∫ δui + ρ 0bi − ρ0 ui dΩ0 + ∑ ∫ 0 δui n0j Pji − ti0 dΓ0 + ∫ δui nj0Pji dΓ0 = 0 (4.8.10)Γ0 ∂X ji= 1 tiΩ0Γint()Since the above holds for all δu∈U0 , it follows by the density theorem given in Section 4.3.2that the momentum equation (B4.5.2) holds on Ω 0 , the traction boundary conditions (B4.5.8)0hold on Γt0 , and the interior continuity conditions (B4.5.11) hold on Γint.
Thus the weak formimplies the momentum equation, the traction boundary conditions, and the interior continuityconditions.4.9 FINITE ELEMENT SEMIDISCRETIZATION4.9.1.Discrete Equations. We consider a Lagrangian mesh with the same properties asdescribed in Section 4.4.1. The finite element approximation to the motion is given byxi (X,t ) = xiI ( t )N I ( X)(4.9.5)where NI (X) are the shape functions; as in the updated Lagrangian formulation, the shpaefunctions are functions of the material (Lagrangian) coordinates.
The trial displacement field isgiven byui ( X, t ) = uiI ( t) N I (X )oru( X, t ) = u I ( t) NI (X )(4.9.1)The test functions, or variations, are not functions of time, soδui (X) = δuiI N I ( X)δu(X ) = δu I N I (X)or(4.9.2)As before, we will use indicial notation where all repeated indices are summed; upper case indicespertain to nodes and are summed over all relevant nodes and lower case indices pertain tocomponents and are summed over the number of dimensions.Taking material time derivatives of (4.9.1) gives the velocity and accelerationu˙ i ( X, t ) = ˙uiI (t ) N I (X)(4.9.3)˙u˙ i ( X, t) = ˙u˙ iI ( t) N I (X )(4.9.4)The deformation gradient is then given byFij =∂x i ∂NI=x∂X j ∂X j iI(4.9.6)It is sometimes convenient to write the above asFij = B 0jI uiI where B 0jI =∂NI∂X jso F = xB 0T4-65(4.9.7)T.
Belytschko, Lagrangian Meshes, December 16, 1998δFij =∂NI∂N IδxiI =δu∂X j∂X j iIso δF = δuB0T(4.9.8)where we have used δxiI = δ ( XiI + uiI ) = δuiI . Nodal forces will now be developed for each of thevirtual weak terms.Internal nodal forces. The internal nodal forces are defined in terms of the internal virtual workusingδW int = δuiI fiIint =∫ δFij Pji dΩ0 = δuiI ∫Ω0Ω0∂NIP dΩ∂X j ji 0(4.9.9)where Eq.
(4.9.7) has been used in the last step. Then the arbitrariness of δuiI yieldsfiIint =∂N∫ ∂XIj PjidΩ 0fiIint = ∫ BjI0 Pji dΩ0 ororΩ0Ω0fiIint = ∫ BjI0 Pji dΩ0(4.9.10)Ω0which is identical to (4.7.6), the expression developed by transformation.External Nodal Forces. The external nodal forces are defined by equating the virtual externalwork (B4.6.3) to the virtual work of the external nodal forces:δWext=δuiI fiIext= ∫ δuiρ 0bi dΩ 0 +Ω0∫δuiti0 dΓ0Γt0i0= δuiI ∫ NI ρ0 bidΩ0 + ∫ N I ti dΓ0 (4.9.12a)Ω0Γt0iThis givesfiIext =∫ NI ρ0bi dΩ0 + ∫ NI ti dΓ00Ω0(4.9.12b)Γt0Mass Matrix: Using the inertial force (B4.6.4) and defining an equivalent nodal force givesδ M =δuiI fiIinert =∫ δui ρ 0˙u˙ idΩ 0(4.9.13)Ω0Substituting Eq. (4.9.2) and (4.9.4) in the right hand side of the above givesδuiI fiIinert = δuiI∫ ρ0 NI NJ dΩ 0˙u˙ jJ = δuiI MijIJ˙u˙ jJΩ0Since the above holds for arbitrary δu and ˙u˙ , it follows that the mass matrix is given by4-66(4.9.14)T.
Belytschko, Lagrangian Meshes, December 16, 1998MijIJ =δ ij ∫ ρ0 NI NJ dΩ0(4.9.15)Ω0Comparing this mass matrix M to that used for the updated Lagrangian formulation, Eq. (4.4.51),we see that they are identical, which is expected since we transformed the mass to the referenceconfiguration to highlight its time invariance for a Lagrangian mesh.Substituting the above expressions into the weak form, (B4.6.1), we haveδuiI ( fiIint − fiIext + MijIJ˙u˙ jJ ) = 0∀I, i∉Γui(4.9.16)Since the above for arbitrary values of all nodal displacement components that are not constrainedby displacement boundary conditions, it follows thatMijIJ˙u˙ jJ + fiIint = fiIext∀I,i ∉Γu i(4.9.17)The above equations are identical to the governing equations for the updated Lagrangianformulation, as given in Box 5.5. The nodal forces in the updated and total Lagrangianformulations are expressed in terms of different variables and integrated over different domains,but from a fundamental viewpoint the updated Lagrangian formulation and the total Lagrangianformulation are identical.
The numerical values for the nodal forces obtained by either formulationare also identical. Each of these formulations can be advantageous for certain constitutiveequations or loadings by reducing the number of transformations which are needed.4.9.2.Implementation. The procedure for the computation of the internal nodal forces isgiven in Box 4.7. In the procedure shown, the nodal forces are evaluated by numericalquadrature.4-67T. Belytschko, Lagrangian Meshes, December 16, 1998Box 4.7Discrete Equations and Internal Force Computation inTotal Lagrangian FormulationEquations of Motion (discrete momentum equation)MijIJv˙ jJ + f iIint = fiIext for ( I,i ) ∉ΓviInternal Nodal ForcesfiIint = ∫ B 0 Ij Pji dΩ0 =Ω0∫Ω0∂N IP dΩ or∂X j ji 0(B4.8.1)(f intI )T= ∫ B0TI PdΩ0(B4.8.2)∫ N I ρ0 bdΩ 0 + ∫ N I ei ⋅ tdΓ0 (B4.8.3)Ω0Tf intI = ∫ B0 I { S}dΩ 0 in Voigt notationΩ0External Nodal ForcesfiIext = ∫ N I ρ0 bidΩ 0 +Ω0∫ N I ti dΓ00or f Iext =Γt0i0Ω0Γt0iMass Matrix (total Lagrangian)MijIJ =δ ij ∫ ρ0 NI NJ dΩ0 = δij ∫ ρ0 N I N J Jξ0d∆Ω0(B4.8.4)∆˜ = I ρ N N dΩM IJ = IMIJ∫ 0 I J 0(B4.8.5)Ω0Internal nodal force computation for element1.
f int = 02. for all quadrature points ξQ[ ] [ ( )]i. compute B Ij0 = ∂N I ξ Q ∂X j for all I∂NIu∂X j iIIiii. F = I + H, J = det( F )1iv. E = (H + H T + H T H )2˙ =∆E ∆t, ˙F =∆F ∆t , D = sym F˙ F-1v. if needed, compute Eii. H = ∑ B 0 I uI ;Hij =((B4.7.6)(B4.7.7))vi.
compute the PK2 stress S or Cauchy stress σ by constitutive equationvii. P = SF T or P = JF−1σintT0viii. f intI ← fI +BOI PJξ wQ for all nodes Iend loopwQ are quadrature weightsUsually the shape functions are expressed in terms of element coordinates ξ , such as the areacoordinates in triangular elements or reference coordinates in isoparametric elements. Thederivatives with respect to the material coordinates are then found by4-68T. Belytschko, Lagrangian Meshes, December 16, 1998( )0N I ,X ≡ BI0 = N,ξ X−1,ξ = N,ξ Fξ−1(4.9.18)where Fξ0 is the Jacobian between the material and intrinsic coordinates. As shown in Box 4.7,the Green strain tensor is usually not computed in terms of the deformation gradient F, because theresulting computation is susceptible to round-off errors for small strains.
Therefore the procedureshown in Eqs. (B4.7.6-7) is used. The total Lagrangian formulation can easily be adapted toconstitutive equations expressed in terms of the Cauchy stress: it is only necessary to introduce thetransormations shown in steps 2.vi-vii.Voigt Form. It is of little use to write the nodal forces in terms of P using Voigt notation since P isnot symmetric. Therefore, we will write the Voigt form in terms of the PK2 stress S . Using thetransformation P = S⋅F T , the expression for the internal nodal forces becomesf jIint =∂N I∫ ∂X FS dΩ 0jk ikΩ0or(f )int TIi=∂NI∫ ∂X SF dΩT(4.9.19)0Ω0We define a B0 matrix by ∂N0BikIj= sym(i, k ) I Fjk ∂Xi(4.9.20)Note that the above specializes to the updated form (4.5.21) where the current configuration andthe reference configuration coincide, so that Fij →δ ij . The Voigt form of this matrix (seeAppendix A) is0BikIj→ Bab0(i, k) → a by the Voigt kinematic rule( I, j) → b by the rectangular to column matrix rule(4.9.21)Similarly, Sik is converted to Sb by the kinetic Voigt rule.
Thenf aint =∫ (B ) S dΩ0 TabbΩ00or f = ∫ B0T {S}dΩ0 or f I = ∫ BT0 I {S}dΩ0Ω0(4.9.22)Ω0The construction of the B0 matrix hinges on the correspondence between the index a and thei and k indices given in Table A.?. Using this correspondence for a two-dimensional element, weobtain:4-69T. Belytschko, Lagrangian Meshes, December 16, 1998B0ikIj → B0 aIj[ ]I∂N I∂N ∂xFj1 = I j∂X∂X ∂X∂N∂N ∂x= I Fj2 = I jI∂Y∂Y ∂Y∂N∂N∂N ∂x ∂N ∂x= I Fj2 + I Fj1 = I j + I jI∂X∂Y∂X ∂Y ∂Y ∂Xi = 1, k = 1 → a =1 B0aj =i = 2, k = 2 → a = 2[ ]B0 aj[ ]i =1, k = 2 → a = 3 B0 aj(4.9.23)The B0 I matrix is then written out by letting j=1 and 2 correspond to columns 1 and 2 of thematrix, respectively:B0 I = ∂N I ∂X∂N I ∂x∂X ∂X∂NI ∂x∂Y ∂Y∂x ∂NI ∂x+∂Y ∂Y ∂X∂N I ∂y∂X ∂X∂NI ∂y∂Y ∂Y∂NI ∂y ∂NI ∂y+∂X ∂Y ∂Y ∂X (4.9.24)In three dimensions, a similar procedure yieldsB0 I = ∂NI ∂Y ∂NI ∂X∂NI ∂X∂N I ∂x∂X ∂X∂NI ∂x∂Y ∂Y∂NI ∂x∂Z ∂Z∂x ∂NI+∂Z ∂Z∂x ∂NI+∂Z ∂Z∂x ∂NI+∂Y ∂Y∂x∂Y∂x∂X∂x∂X∂N I∂Y∂NI∂X∂NI∂X∂N I ∂y∂X ∂X∂NI ∂y∂Y ∂Y∂NI ∂y∂Z ∂Z∂y ∂NI+∂Z ∂Z∂y ∂NI+∂Z ∂Z∂y ∂NI+∂Y ∂Y∂y∂Y∂y∂X∂y∂X∂N I ∂z∂X ∂X∂NI ∂z∂Y ∂Y∂NI ∂z∂Z ∂Z∂N I ∂z ∂N I+∂Y ∂Z ∂Z∂N I ∂z ∂NI+∂X ∂Z ∂Z∂NI ∂z ∂NI+∂X ∂Y ∂Y∂z ∂Y ∂z ∂X ∂z ∂X (4.9.25)Many writers construct the B0 I matrix through a sequence of multiplications by Boolean matrices.The procedure shown here can easily be coded and is much faster.˙ to the node velocities byIt can be easily shown that B0 relates the rate of Green strain E{E˙ } = B0I v I = B0˙d(4.9.27)The reader should be cautioned about one characteristic of the B0 matrix: although it carries asubscript nought, the matrix B0 is not time invariant.