Belytschko T. - Introduction (779635), страница 76

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 76)

As a result, we may select either formulationdepending on which is more convenient computationally for the problem at hand.d 2φFinally, it can be noted that when linear elements are used, 2 = 0 , so all terms involvingdxsecond derivatitves vanish. Consequently, Eq. (7.14.31) may be written as:0 =∑e∫Ωe wu dφ +v ∗ dw dφ  dxdxdx dx (7.14.32)where, υ* ,which may be thought of as a sum of two viscosities will be defined below.7.14.4 Parameter Determination and Further AnalysisTo begin, we wish to shed additional light on the physical interpretation of Eq. (7.14.32).As a result, we make the following definitions:υ * = υ + υ = total viscosityandυ = αu∆xsign(u), α ≥ 02(7.14.33a)(7.14.33b)W.K.Liu, Chapter 720It then becomes clear that υ may be thought of as an artificial viscosity which must beadded to the “normal” flow viscosity, υ , to ensure stability. That is, without thissuperficial damping which does not correspond to the physics of the problem, ournumerical solution oscillates wildly thereby leading to physically meaningless results.To define these viscosities in terms of various Peclet numbers, consider the followingrelationships:υ ⇐ Pe =u∆x2υ(7.14.34a)υ ⇐ Pe =u∆x2υ(7.14.34b)υ * ⇐ Pe* =u∆x2υ *(7.14.34c)The relationships between Eq.

(7.14.34a-c) can be expressed as:12υ * 2υ2υ11==+=+*Peu∆x u∆x u∆x Pe Pe(7.14.35)In Eq. (7.14.35), if Pe* < 1, then2υ2υ11+> 1 or> 1−u∆x u∆xPePe(7.14.36)and the solution will not be oscillatory. From Eq. (7.14.13), the discrete equation in termsof Pe* can be written as:(Pe * − 1)φ j+1 + 2φ j − (Pe* + 1)φ j−1 = 0(7.14.37)Recall that the discrete solution is oscillatory when φ N = z N , the roots of Eq. (7.14.37)are : 1 + Pe*  Nz = c1 , c2  1 − Pe* N(7.14.38)From Eq. (7.14.21) and Eq. (7.14.38), we can solve Pe* .

That is,u( ∆x) j 1 + Pe*  jυ=e 1 − Pe* (7.14.39)W.K.Liu, Chapter 721The RHS of Eq. (7.14.39) can be expressed in terms of Pe . That is:uu∆x∆x2()1 + Pe*υ2υ = e2 Pe=e=e1 − Pe*(7.14.40a)1 − e2 Pe = −Pe* (1+ e 2 Pe )(7.14.40b)orTherefore,e 2P e − 1 e Pe − e − PePe = 2P e=e+ 1 e Pe + e − Pe(7.14.41a)Pe* = tanh(Pe )(7.14.41b)*orSubstituting Eq. (7.14.41b) into Eq. (7.14.35) yields:111=+tanh( Pe ) Pe Pe(7.14.42a)11= coth( Pe ) −PePe(7.14.42b)orTogether, Eq. (7.14.34b) and Eq. (7.14.42b) may be combined to give:u∆x 1Pe ==  coth(Pe ) − 2υ Pe −1(7.14.43)Using Eq.

(7.14.43), we can also express υ in terms of Pe :υ =11 u∆xu∆x  coth( Pe ) −  = αsign(u)2Pe 2so therefore,(7.14.44a)W.K.Liu, Chapter 7υ =2211 u∆xu∆x  coth( Pe ) −  = αsign(u)2Pe 2(7.14.44b)Finally, it becomes apparent that we may define the parameter α as:α = sign( u) [coth( Pe ) −1]Pe(7.14.45)Note that when α = 0 , Eq. (7.14.29) is simply like a central difference method and whenα = 1, it is a full upwind Petrov-Galerkin formulation.7.14.5 Streamline-Upwind/Petrov-Galerkin Formulation for Multiple DimensionsThe advection-diffusion equation in multiple dimension is:u ⋅∇φ − v∇ 2 φ = 0in Ω(7.14.45a)where the boundary conditions are:φ=gon ∂Ω g(7.14.45b)v ∇φ ⋅ n = 0on ∂Ωt(7.14.45c)The weak form of the advection-diffusion equation for a streamline-upwind/PetrovGalerkin formulation of Eq. (7.14.45a) is obtained by multiplying Eq.

(7.14.45a) by thetest function w and integrating over the domain Ω∫ w( u⋅∇φ − v∇ φ) = 02Similar to the one-dimensional formulation presented above, let the Petrov-Galerkin testfunctions be defined by˜≡w →ww+τu⋅∇w{123Galerkintestfunctiondiscontinuoustest function(7.14.45d)thus giving0=∫Ω(w + τu ⋅ ∇w)(u ⋅∇φ − υ∇ 2 φ)dΩwhere the stabilization parameter is now given by(7.14.46a)W.K.Liu, Chapter 7τ=α23hand h ≡ ∆x2u(7.14.45b)where α is given by Eq.(7.14.33c). For a time dependent problem, the stabilizationparameter can be set by:τ=α∆t2(7.14.45c)˜ ∉U 0 as in the one dimensional case. Applying integration by partsNote that w ∈U 0 and wand the divergence theorem, the weak form of Eq. (7.14.46a) can be shown to be:0 = ∫ wu ⋅∇φ dΩ+ ∫ υ∇w ⋅ ∇ φdΩΩ44444Ω 4444414243Galerkin termsNeNe+ ∑ ∫ e τ (u ⋅ ∇ w)(u ⋅∇φ) dΩ− ∑ ∫ e τ (u ⋅ ∇w)υ∇ 2 φdΩΩΩ=1 4444444441e=14444444444424e43(7.14.46)streamlineupwind stabilization termsAs can be seen from the above equation, the Petrov-Galerkin terms are simply the sum ofthe standard Galerkin terms plus the streamline upwind stabilization terms.

Namely,Petrov / Galerkin = Galerkin + Streamline Upwind StabilizationThe third term of Eq. (7.14.46) can be rewritten as:∫Ωe=τ (u ⋅∇w)(u ⋅∇φ) dΩ∫ τ [wΩe,1=∫ τ [w=∫whereΩeΩe,1w,2 φ,1  u1 w,3 u2  [u1 u2 u3 ] φ ,2  dΩφ  u 3 ,3 w,2 u1u1 u1u2w,3 u2 u1 u2u2 u u u u3 13 2]∇w ⋅ (υ˜ ∇φ) dΩ]u1 u3   φ ,1 u2 u3   φ,2  dΩu3u3   φ,3 (7.14.47a)W.K.Liu, Chapter 724 u1u1 u1u2υ˜ = τ u2 u1 u2u2 u u u u3 13 2u1u3 u2u3 u3u3 (7.14.47b)Substituting Eq.

(7.14.47a) into Eq. (7.14.46) and ignoring the last second order term, Eq.(7.14.46) becomes:artifical}viscosity0 = ∫Ω w u ⋅ ∇φ dΩ + ∑ ∫Ω e ∇w(υ1 +υ˜)∇φ dΩe =11444442444443Galerkin TermNe(7.14.48)The artificial viscosity acts as a stabilization term which eliminates the oscillations resultingfrom a standard Galerkin formulation.7.14.6 An Alternative Derivation of the Multiple Dimensional Streamline-Upwind/PetrovGalerkin FormulationParalleling section 7.14.3, this section outlines an alternative derivation of themultidimensional streamline-upwind/Petrov-Galerkine formulation presented above. Tothis end, we begin by starting with the multiple dimensional advection-diffusion equationas before. The equation is then multiplied by a test function w˜ and integrated over thedomain Ω (following the traditional weak form development) thereby yielding∫˜ ( u⋅∇φ − v∇2 φ) = 0wExamining the second term on the left hand side in more detail, our goal is to remove onederivative from the trial function and place it on the test function, w˜ , thereby relaxing trialfunction continuity requirements.

To this end, we integrate by parts in the familiar manneras follows:˜ v∇ 2φdΩ = ∫ ∇( w˜ v ∇φ) dΩ − ∫ ∇w˜ v∇φdΩI ≡ ∫Ω wΩΩ˜ gives:Applying the divergence theorem and the substituting in the definition of w˜ v∇φ Γ − ∫ ∇ w˜ v∇φdΩ= [w +τu⋅∇w]v∇φI=wΩ= τu ⋅∇ wv∇φΓΓ˜ v∇φdΩ− ∫ ∇wΩ˜ v∇φdΩ− ∫ ∇wΩwhere the last line has been obtained by using the fact that w ∈U 0 .

Combining theabove results with the advection term yields the following alternative weak form:∫ [w˜ u⋅∇w +∇ w˜ v∇φ ]dΩ− τu ⋅∇wv∇φΩΓ=0Now it is apparent that removal of a trial function derivative gives rise to a boundaryintegral term which was not present in the streamline-upwind/Petrov-Galerkin formulationW.K.Liu, Chapter 725presented earlier.

In the particular case when w˜ is defined as in Eq. (7.14.28), it isstraightforward to show that this alternative formulation yields the same results as theformulation presented in the previous section.Example Petrov-Galerkin formulation of the ALE momentum equation for a 1D 2-nodelinear displacement element may be done as follows.The test function is chosen to be:h dN INI = N I + αsign(c)2 dxwhere α is the viscosity constant, c is the convective velocity, h is the element size, andNI is the shape function of the usual Galerkin form. If constant density is assumed, weobtain the Petrov-Galerkin form of the mass matrix in the ALE formulation, as: 1 − α (K − K ) 1 − α K 21126 2 M = ρh  31 α1 α + (K2 − K1 )+ K 6 23 2 1where 1 sign(c )122K1= 1 c1  + 1 sign(c2 ) [sign(c1 ) − sign(c2 )] c1 − c2 22sign(c1 )K 2 =  sign(c ) − sign(c )  c1  + sign(c )12 ]2[ c1 − c2 if (c1 = c2 )if (c1 ≠ c2 )if (c1 = c2 )if (c1 ≠ c2 )and c1 and c2 are the convective velocities at the 2 nodes.Similar to that for convective-advective equation, we can rewrite the mass matrix as: 1 1 − α (K − K ) − α K 21122 M = ρh  31 61  + ρh  αα + (K2 − K1 ) + K1  426 4332 441414424442443Galerkintermsstreamline upwind stabilizationtermsnamely, the Petrov-Galerkin terms can be split into the sum of Galerkin terms and thestreamline upwind stabilization terms.7.3 Weak Form and Petrov-Galerkin Finite Element Discretization of theALE Continuity and Momentum Equations:As can be seen from Box 7.1, the only difference between the updated Lagrangianequations and the updated ALE equations is the interpretation of the material timederivative.

Hence the weak form, and subsequently the Galerkin finite elementW.K.Liu, Chapter 726formulations, are identical to those derived in Chapter 4. However, note that the spatialdomain now depends on how the mesh motion is updated, which is one of the keyingredients in the updated ALE formulation.The variational equations corresponding to the conservation equations of Box 7.1 areobtained by multiplying by the test functions, δρ and δvi , integrating over the spatialdomain Ω and employing the divergence theorem to embed the traction force vector t onthe boundary Γ t .

Following the same procedures as in Chapter 4, we can achieve thefollowing weak forms:Continuity Equation:∫Ωδρρ˙dΩ = ∫Ω δρρ ,t [χ] dΩ+∫ δρc ρiΩ,idΩ = − ∫Ω δρρvi,i dΩ(7.16.2a)Momentum Equation:∫Ωδviρ v˙i dΩ= ∫Ω δvi ρvi,t [χ ]dΩ + ∫Ω δv iρc j vi, j dΩ= − ∫ δvi, j σ ij dΩ+Ω(7.16.2b)∫ δv ρb dΩ + ∫ΩiixΓtδvit idΓIt is noted that because convective terms (ρ,i ci and vi, j c j ) appeared in the continuity andmomentum equations, a Galerkin finite element formulation will give rise to numericaldifficulties.

Therefore, in this section, the Petrov-Galerkin formulation will be employed toalleviate some of these difficulties. In a Petrov Galerkin finite element discretization, thecurrent domain Ω is subdivided into elements. However, different sets of shape functions,N and Nρ for the trial functions, and N and Nρ for the test functions, will be used tointerpolate the velocity and density, respectively. If N = N , the Galerkin ALE formulationswill be obtained.

The choice of N and Nρ to eliminate numerical oscillations will bedescribed in section ????.The finite element matrix equations corresponding to Eq.(7.16.2a,b) are :Continuity equation:M ρρ,t[χ ] + Lρ ρ + Kρ ρ = 0(7.16.3a)where M ρ , Lρ , K ρ are generalized mass, convective, and stiffness matrices,respectively, for density under a reference description such that:M ρ = [MIJρ ] = ∫ N Iρ N ρJ dΩ(7.16.3b)Lρ = [LρIJ ] = ∫ NIρ ci NJρ,i dΩ(7.16.3c)ΩΩW.K.Liu, Chapter 727K ρ = [KIJρ ] = ∫ NIρ vi,i N ρJ dΩ(7.16.3d)ΩMomentum Equation:Mv ,t[ χ ] + Lv + f int = f ext(7.16.4a)where M and L are generalized mass and convective matrices, respectively, for velocityunder a reference description; while f int and f ext are the internal and external force vectorsrespectively, such that:M = [M IJ ] = ∫ ρNI N J dΩ(7.16.4b)L = [LIJ ] = ∫ ρNI ci N J ,i dΩ(7.16.4c)f int = [ f iIint ] = ∫ NI , j σ ij dΩ(7.16.4d)ΩΩΩf ext = [ f iIext ] = ∫ ρNI bi dΩ + ∫ t N I ti dΓΩΓ(7.16.4e)REMARK 7.16.1.4 The nonself-adjoint and nonlinear convective terms, Lρ and L ,which appear in Eqs.(7.16.3c) and (7.16.4c) and characterize the ALE method, willinevitably pose difficulties.REMARK 7.16.1.5 All the matrices and vectors defined in Eqs(7.16.3)-(7.16.6) areintegrated over the spatial domain Ω .

Характеристики

Список файлов книги