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

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

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

A closed form solution for thediscrete steady-state advection diffusion equation will be obtained. It will be shown that thesolution is oscillatory when a parameter of the mesh, known as the Peclet number, exceedsa critical value. Next, a Petrov Galerkin method will be developed which eliminates theseoscillations.Consider the linear advection-diffusion equation:∂φ+ u ⋅ ∇φ − υ ∇ 2 φ = 0∂t(7.14.1)where φ is the dependent variable, υ is the kinematic viscosity, and u is a given velocity∂φfield. For the steady state case,= 0. So, the steady state equation is:∂tu ⋅∇φ − υ ∇ 2 φ = 0(7.14.2)For the study of special numerical instabilities, we restrict Eq.

(7.14.2) to one dimensionso that:udφd 2φ=υ 2dxdx(7.14.3)Equation (7.14.3) with boundary conditions:φ (0) = 0 and φ (L) = 1is a two-point boundary value problem on the domain 0 ≤ x ≤ L .It is easy to verify that the exact solution to Eqs. (7.14.3) and (7.14.4) is:(7.14.4)W.K.Liu, Chapter 7φ (x) =141 − e ux/ υ1 − euL /υ(7.14.5)7.14.1 The Galerkin Finite Element Approximation of the Advection-Diffusion EquationLetting the test function be w(x), multiplying Eq.

(7.14.3) by w, and integrating overthe domain gives dφd 2φ w u− υ 2  dx = 0dx Ω  dx∫(7.14.6)Integrating by parts and making use of the divergence theorem, the weak form of the onedimensional advection-diffusion equation, Eq. (7.14.3), is∫Ωwudφdx +dx∫ υw φΩdx = 0,x ,x(7.14.7)with w ∈U 0 . The domain (0, L) is then divided into equally sized linear finite elements,Ω e , on which the finite element approximation is given by:(∫Ωeu N a Nb,x dx)φ b + (∫Ωeυ Na,x N b,x dx)φ b = 0a, b = 1, 2(7.14.8)where Na and Nb are the linear finite elements shape functions.

This can be written inmatrix component form asNab φb + Kab φ b = 0a, b = 1, 2(7.14.9a)where the convective matrix is given as:Nab =∫x e+1u N a Nb,x dxxe(7.14.9b)and the diffusion matrix is:Kab =∫x e+ 1xeυ N a,x Nb,x dx(7.14.9c)It is a simple exercise to show that, when using linear finite element shape functions,N=u  −1 12  −1 1K=υ∆x 1 −1 −1 1 After assembly, the equation for the jth node is :(7.14.10)W.K.Liu, Chapter 715 φ j +1 − φ j −1  φ j+1 − 2φ j + φ j −1 u− υ=02∆x ∆x 2(7.14.11a)which is exactly the central difference approximation.

It is convenient to normalize theabove equation, so that:() ()u∆xφ j+1 − φ j−1 − φ j+1 − 2φ j + φ j −1 = 02υ(7.14.11b)The Peclet number, Pe , may then be defined as:Pe =u∆x2υ(7.14.12)In terms of the Peclet number, Eq. (7.14.11b) then becomes:() ()Pe φ j +1 − φ j −1 − φ j+1 − 2φ j + φ j−1 = 0( Pe − 1) φ j +1 + 2φ j − ( Pe + 1) φ j −1 = 0(7.14.13a)(7.14.13b)Ignoring the boundary conditions, Eq. (7.14.13) can be put into the standard matrixnotation by expanding the jth term in Eq.

(7.14.13): −1 0 1−1 0 1Pe −1 0convective term M    M  M  φ j −1   −1 2 −1  φ j−1   0   φj  +   φ j  = 0 −1 2 −11 φ j +1  −1 2 −1  φ j+1   0     M    M  M diffusion term(7.14.14)The solution to the discrete finite difference Eq. (7.14.13) can be obtained by assuming:φ (x j ) ≡ φ j = eax j= e a( j ∆x ) = e ( a∆x ) j ≡ z j(7.14.15)Where z = e a∆x and a is an unknown coefficient to be determined. By the definition in Eq.(7.14.15), the j+1th and j-1th terms of φ are:φ j +1 = e a( j+1) ∆x = e a j∆x e a∆x = z j +1(7.14.16a)φ j −1 = e a( j−1)∆x = e a j∆x e −a ∆x = z j−1(7.14.16b)Substituting Eqs. (7.14.16) into Eq.

(7.14.13) yields:W.K.Liu, Chapter 716( Pe − 1) z j +1 + 2z j − ( Pe + 1) z j −1 = 0(7.14.17)Assuming that z j−1 ≠ 0 and dividing the above equation by z j−1 , Eq. (7.14.17) becomes:( Pe − 1) z2 + 2z − (Pe + 1) = 0(7.14.18)The roots for Eq. (7.14.18) are:z = 1 or z =1 + Pe1 − Pe(7.14.19)Recalling that φ j = z j , the solution to Eq. (7.14.13) takes the form: 1 + Pe  jφ j = c1 + c2  1 − Pe (7.14.20)where c1 and c2 are coefficients to be determined from the boundary conditions. Since theexact solution to Eq. (7.14.3) is given by Eq.

(7.14.5), the exact solution of φ , evaluatedat x = x j , has the form of:[]uj∆x1ux j / υφ (x j ) == c1 + c2 e υuL/ υ 1− e1− e(7.14.21)Comparing the finite difference solution Eq. (7.14.20) with the exact solution, Eq.(7.14.21), it can be concluded that:(i) If the Peclet number is less than one, i.e., Pe < 1, then the discrete solution willhave a solution similar to the exponential solution as given in the exact solution since 1 + Pe  j > 0. 1 − Pe (ii) If the Peclet number is greater than one, i.e., Pe > 1, then the discrete solutionbecomes: 1 + Pe  j = (−m) j with m > 0 1 − Pe Hence nodal oscillations occur because φ j is positive or negative depending on whether jis even or odd, respectively.

To illustrate these nodal oscillations, we consider the oneW.K.Liu, Chapter 717dimensional advection-diffusion equation as given in Eq.(7.14.3) with boundaryconditions (7.14.4). The plots below compare the exact solution with finite elementsolutions for the cases of both no upwinding and full upwinding. In all cases, 80 elementswere used with an element Peclet number of 300.7.14.2 Ramification of Nodal Oscillation by the Petrov-Galerkin FormulationRecall the weak form of Eq.

(7.14.3):∫Ωdφd2 φw (u− υ 2 )dx = 0dxdx(7.14.27)The Petrov-Galerkin formulation for Eq. (7.14.3) is obtained by replacing the test functionw by w˜ , where w˜ is defined as:˜ ≡w→w+w{Galerkintestfunction∆x dwα(signu)24dx142443discontinuoustestfunction(7.14.28)Replacing w by w˜ , Eq. (7.14.27) becomes:∫Ωw˜ (udφd2 φ− υ 2 )dx = 0dxdx(7.14.29)Note that w ∈U 0 and w˜ ∉U 0 . The parameter α is to be determined so as to eliminateoscillations for Pe > 1 and hopefully get accurate solutions; in one dimension, it is possibleto select α so as to obtain exact values of the solution at the nodes. Substituting thedefinition of w˜ , Eq.

(7.14.28), into Eq. (7.14.29) yields:0 = ∫Ω w˜ (udφd 2φdφd 2φ− υ 2 )dx = ∫Ω w(u− υ 2 )dxdxdxdxdx44144442443GalerkinTermNe∆x dwdφd 2φ+ ∑ ∫Ωe α(signu)(u− υ 2 )dx2 4dxe=11444444244dx444dx443Upwind Petrov − Galerkin Term(7.14.30)After integrating by parts (and using w (0)= w (L)=0 be construction), Eq. (7.14.30)becomes:0=∫ wu dx dx + ∫ υ dx dx dxdφΩdw dφΩW.K.Liu, Chapter 7Ne+∑∫e =1Ωe18α∆xdw dφu(sign u)dx −2dx dxNe∑∫e=1Ωeυα∆xdw d 2 φ(signu)dx2dx dx 2(7.14.31)The above equation is known as upwinding Petrov-Galerkin formulation(Brooks &Hughes, 1978).

It is noted that in this formulation the second derivative of φ is required.Further, the free parameter α is determined in the following section after the presentationof an alternative formulation which only requires the first derivative of φ .7.14.3 An Alternative Derivation of the Upwind/Petrov-Galerkin FormulationThis section outlines an alternative derivation of the upwind formulation. Motivated by adesire to necessitate low order continuity restrictions on the trial functions, we begin bystarting with the one dimensional advection-diffusion equation as before. The equation isthen multiplied by a test function w˜ and integrated over the domain Ω (following thetraditional weak form development) thereby yieldingdφd 2φ˜w(u−υ)dx = 0∫Ωdxdx2Examining 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:I ≡ φ ∫Ω˜d2φ˜˜ υ dφ ) , x dx − ∫ d w υ dφ dxw υ 2 dx = ∫Ω ( wΩdxdxdx dx˜ gives:Applying the divergence theorem and the substituting in the definition of wL˜ dφ˜ dφdφ Ldw∆x dwdw dφ− ∫Ωυdx = w + α(sign u) υ− ∫Ωυdx dx 0dx odx dx2 dxdx dx˜ dφ∆x dwdφ Ldw=α(sign u)υ− ∫Ωυdx( since w( o) = 0 and w( L) = 0)2 dxdx 0dx dxCombining the above results with the advection term yields the following alternative weakform:˜υI=w˜ dφ ∆x dwdφ L ˜ dφ dwwu+υdx−α(signu)υ=0∫Ω  dx dx dx 2 dxdx 0Now it is apparent that removal of a trial function derivative gives rise to a boundaryintegral term which was not present in the Petrov-Galerkin formulation presented earlier.In the particular case when w˜ is defined as in Eq. (7.14.28), it is straightforward to showthat this alternative formulation yields the same results as the formulation presented in theprevious section.

To see this, substitute the explicit expression for w˜ into the equationgivingW.K.Liu, Chapter 719∫ dw∆x d 2w dφ w + α ∆x dw (sign u)  u dφ dx ++α(sign u)  vdx∫2Ω2 dx dx2 dx dx dx−α∆x dwdφ(sign u)υ2 dxdxΩL=00Using integrating by parts on the fourth term gives∫Ω ∆x d 2w dφ∆xdw dφ L∆xdw d 2φα 2 dx 2 (sign u)  v dx dx =α 2 (sign u)v dx dx 0 − ∫Ω α 2 (sign u)v dx dx 2 dxwhich, upon rearrangement of the terms, yields the expression resulting from thepreviously presented Petrov-Galerkin formulation:∫ΩNe−∑∫Ωeα∆xdw dφ(sign u)udx2dx dx∆xdw d 2φ∆xdw dφα(sign u)vdx +α(sign u) v22dx dx2dx dx∫Ωee =1αNedφdw dφdx + ∫ vdx + ∑Ωdxdx dxe =1wuL−0∆xdw dφ L(sign u)v=02dx dx 0Following the cancelation of the last two terms, this equation is identical to that of theprevious formulation, Eq.(7.14.31).

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

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