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

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

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

An Euleriandescription used in the x1 and x2 directions (i.e. x1 = χ 1 and x2 = χ 2 ). The free surface isdefined by one spatial coordinate which is a continuous and differentiable function of theother two spatial coordinates and time. In this case the Lagrange Euler matrix has only onenon-zero term, α 33 (usually equal to 1), and the only non-trivial equation in (7.11.3) is:∂x3s∂x∂x∂x+ v1 3s + v2 3s − v3 = (α 33 − 1) v3 3s∂t χ∂χ1∂χ 2∂χ 3(7.11.7)The above equation is easily recognized as the kinematics equation of the surface and maybe written as:∂x3s+ vi ni Ns = a(x1 , x2 ,x 3s ,t)∂t χ(7.11.8)where the components of n is the unit normal pointing out from the surface.components of the normal vector are given by:1  ∂x 3s ∂x3s,,−1Ns  ∂χ1 ∂χ 2The(7.11.9a)with Ns given by:1/2  ∂x  2  ∂x  2 Ns = 1 +  3s  +  3s   ∂χ 2    ∂χ 1 1/2  ∂x  2  ∂x  2 = 1+  3s  +  3s   ∂x2    ∂x1 (7.11.9b)where a(x1 , x2 ,x 3s ,t) is the so-called accumulation rate function expressing the gain or lossof mass under the free surface (Hutter and Vulliet, 1985).

It can be seen by comparing Eq.(7.11.7) and Eq. (7.11.8) that the accumulation rate function is:a(x1 , x2 ,x 3s ,t) = (α 33 − 1)v3∂x3s∂x= w3 3s∂χ 3∂χ 2(7.11.10)The free surface is a material surface; along the free surface the accumulation rate must bezero, and consequently α 33 has to be taken equal to one. This can also be deduced bynoticing that no particles can cross the free surface, so w3 must be zero. Although Eq.(7.11.3) can be applied to problems where x1 and/or x2 are not Eulerian by prescribingnon zero α 's in these directions, controlling the element shapes by adjusting the α 's isvery difficult.W.K.Liu, Chapter 751Control of the mesh by Eq.

(7.11.1) has some disadvantages; for instance, while vˆhas a clear physical interpretation (i.e. the mesh velocity), w is much more difficult tovisualize (except in the direction perpendicular to material surfaces, where it is identicallyzero) and therefore it is very difficult to maintain regular shaped elements inside the fluiddomain by just prescribing the α 's. Because of this important drawback the mixedformulation introduced by Huerta and Liu, called deformation gradient method, wasdeveloped, it is discussed next.7.11.4 Deformation Gradient FormulationsNoticing the restrictions of the α 's scheme, a mixed formulation is developed for theresolution of Eq. (7.2.19). One of the goals of the ALE method is the accurate tracking ofthe moving boundaries which are usually material surfaces.

Hence, along these surfaceswe enforce w ⋅ n = 0 where n is the exterior normal. The other goal of the ALE techniqueis to avoid element entanglement and this is better achieved, once the boundaries areknown, by prescribing the mesh displacements independently (through the potentialequations, for instance) or velocities, because both dˆ and vˆ govern directly the elementshape. Therefore, one can prescribe w ⋅ n = 0 along the domain boundaries while definingthe dˆ ‘s or vˆ ‘s in the interior.The system of differential equations defined in Eq. (7.2.19) has to be solved alongthe moving boundaries.

Notice first that solving for wi in terms of (vi − vˆi ), Eq. (7.2.19)can be rewritten as:c j ≡ v j − vˆ j = F χji wi(7.11.11)Define the Jacobian matrix of the map between the spatial and ALE coordinates byFijχ ≡∂x i∂χ j(7.11.12)Its inverse is:(F )χ −1 Jˆ111= ˆ  − Jˆ21J ˆ J 31− Jˆ12Jˆ22− Jˆ32Jˆ13 − Jˆ23  ≡Jˆ33 JˆijJˆ(7.11.13a)where Jˆij are the cofactors of Fijχ , Jˆ is the Jacobian already defined in Eq.

(7.7.4b) andJˆij are the minors of the Jacobian matrix F χ . Multiplying the inverse Jacobian matrix onijboth sides of Eq. (7.11.11) and substituting Eq. (7.11.13) into Eq. (7.11.11), yields:Jˆijˆ jiˆˆˆˆJ v j − v j = wi or J (v j − v j ) = J wi(7.11.14)W.K.Liu, Chapter 752Dividing Jˆii on both sides of Eqs. (7.11.14), gives: ∂xi− vi i = j ∂t χJˆ jiJˆˆ(v−v)=w= nsd v j − vˆ jjjiJˆ iiJˆiiˆ ji i ≠ jˆJ ii J jj≠=1i∑(7.11.15)When the LHS of Eq. (7.11.15) has been simplified using by the definition of vˆ j , Eq.(7.2.8a). Substituting Eq. (7.11.16) into Eq. (7.11.15) yields:∂xi∂tNSDχ− vi −∑j =1j ≠iv j − vˆ j ˆ j iJˆJ=−wJˆ iiJˆii i(7.11.17)Notice that the cofactor Jˆii appears in the denominator to account for the motion of themesh in the plane perpendicular to χ i because Eqs.

(7.11.17) are verified in the referenceˆ , not in the actual deformed domain Ω .domain ΩExamples for Eq. (7.11.17) in 1D, 2D and 3D are:W.K.Liu, Chapter 7531D∂x1Jˆ− v1 = − ˆ11 w1∂t χJ(7.11.18)where Jˆ11 = 1.2D∂x1v − vˆJˆ− v1 − 2 ˆ11 2 Jˆ 21 = − ˆ11 w1∂t χJJ(7.11.19a)∂x2∂t(7.11.19b)χ− v2 −v1 − vˆ1 ˆ12JˆJ=−wJˆ 22Jˆ 22 2∂x∂xwhere Jˆ11 = 2 and Jˆ22 = 1 .∂χ 2∂χ13D∂x1v − vˆv − vˆJˆ− v1 − 2 ˆ11 2 Jˆ 21 − 3 ˆ11 3 Jˆ31 = − ˆ11 w1∂t χJJJ(7.11.20a)∂x2∂t(7.11.20b)χ− v2 −v1 − vˆ1 ˆ12 v3 − vˆ3 ˆ32JˆJ−J=−wJˆ 22Jˆ 22Jˆ22 2∂x3v − vˆv − vˆ− v3 − 1 ˆ 33 1 Jˆ13 − 2 ˆ33 2 Jˆ 23 = −∂t χJJJˆwJˆ 33 3(7.11.20c)For purposes of simplification and without any loss of generality, assume that themoving boundaries are perpendicular to one coordinate axis in the reference domain. Letthe free surface be perpendicular to χ 3 , the first two equations in Eq.

(7.11.17) are trivialbecause in the direction of χ 1 and χ 2 the mesh velocity is prescribed and therefore themesh motion is known, but the third one must be solved for vˆ3 given w3 , vˆ1 , and vˆ2 , itmay be written explicitly as:Jˆ13Jˆ23Jˆvˆ3 − ˆ 33 (v1 − vˆ1 ) − ˆ33 (v2 − vˆ2 ) − v3 = − ˆ 33 w3JJJ(7.11.21a)W.K.Liu, Chapter 754or ∂x ∂x  v − vˆ ∂x ∂x ∂x3sv − vˆ− 1ˆ 33 1 Jˆ13  3s , 3s  − 2 ˆ 33 2 Jˆ23  3s , 3s  − v3∂t χJ ∂χ 1 ∂χ 2 J ∂χ1 ∂χ 2 w= − ˆ333J ∂x ∂x Jˆ 3s , 3s  ∂χ1 ∂χ 2 where vˆ3 has been substituted by(7.11.22b)∂x3s, Jˆ13 , Jˆ23 , and the Jacobian Jˆ are function of∂t χ∂x3s∂xand 3s : Jˆ33 is not dependent on x3s : and x3s is the free surface equation.

In Eq.∂χ1∂χ 2(7.11.22b) x3s is the unknown function, while vˆ1 , vˆ2 and w3 are known. If vˆ1 = vˆ2 = 0(i.e. the Eulerian description is used in χ 1 and χ 2 ), the kinematic surface equation, Eq.(7.11.8), is again obtained. However, with the mixed formulation vˆ1 and vˆ2 can beprescribed (as a percentage of the wave celerity, for instance) and therefore better numericalresults are obtained than by defining in Eq. (7.11.7) α 11 and α 22 , whose physicalinterpretation is much more obscure.7.11.5 Automatic Mesh GenerationThe Laplacian method for remeshing is based on mapping the new position of thenodes by solutions of the Laplace equation space (I, J) into real space (x, y) throughsolving the Laplace differential equation is the most commonly use one.The determination of the nodes is posed as finding x(I, J) and y(I, J), such that theysatisfy the following equations∂ 2 x ∂ 2x∂2 y ∂2 y2L (x) = 2 + 2 = 0; L (y) = 2 + 2 = 0∂I∂J∂I∂J2in Ω x(7.11.23a)in Γ x(7.11.23b)The boundary conditions in two dimension arex(I, J) = x (I, J); y(I,J ) = y (I, J )Here x(I) and y(J) are the coordinates of nodes I and J when I and J take on integer values.Ω and Γ are the domain and boundary of the remesh region and L2 is the second-orderdifferential operator.Another useful mesh generation scheme is by solving a fourth order differentialequationL4 (x) =∂4 x∂ 4x4;L(y)=∂I 2 ∂J 2∂I 2 ∂J 2(7.11.24)W.K.Liu, Chapter 755Eqs.

(7.11.23) and (7.11.24) can be solved by the finite difference method with a GaussSeidel iteration scheme. Meshes generated by the Laplace equation are distorted near theboundary where a high curvature occurs. However, the fourth-order equation gives abetter mesh shape, because a higher differentiation is employed. An equipotential methodregards the mesh lines as two intersecting sets of equipotentials, with each set satisfyingLaplace’s equation in the interior with adequate boundary condition.7.12 Strong Form, Governing Equations of Slightly Compressible ViscousFlow with Moving Boundary ProblemThe object here is to find the following functions:v(x,t) = velocity, p(x,t) = pressure fields, xˆ (x,t) = mesh displacement(7.12.1c)such that they satisfy the following state and field equations:Continuity Equation [Eq. (7.10.3b)]1 ∂p1 ∂p ∂vi+ ci+=0B ∂t χ B ∂xi ∂xi(7.12.2a)Momentum Equation [Eq.

(7.10.1b)]ρ∂vi∂tχ+ ρc j∂vi ∂σ ij=+ ρgi∂x j ∂x j(7.12.2b)Free Surface Update EquationWe can employ either mesh rezoning equation [Eq. (7.11.3)]:∂xi∂tχ()+ δ jk − α jk vk∂xi− vi = 0∂χ j(7.12.2c)or the mesh updating equation [Eq. (7.11.17)]:∂xi∂tNSDχ− vi −∑j =1j ≠1v j − vˆ j ˆ j iJˆJ=−wJˆ iiJˆii i(7.12.2d)The boundary conditions are as follows. It is required that:vi = vion Γ vi(7.12.3)σ ij n j = tion Γ ti(7.12.4)W.K.Liu, Chapter 756where b and h are the prescribed boundary velocities and tractions, respectively; n is theoutward normal to Γ vi , and Γ vi is the piecewise smooth boundary of the spatial domain, Ωand the decomposition of Γ is given in Chapter 3.Constitutive Equationσ ij = − pδij + 2µDij(7.12.5a)whereDij =1  ∂vi ∂v j +2  ∂x j ∂xi and µ = µ (Dij )(7.12.5b)µ is the dynamic viscosity which is shear rate dependent.

In Table 7.2 several GeneralizedNewtonian models (see e.g. Bird et al, 1977) are presented. The finite element methodpresented here is independent of the particular Generalized Newtonian model chosen.7.13 Weak Form of SlightlyBoundary ProblemCompressibleViscousFlowwithMovingWe denote the spaces of the test function and trial functions by:{δvi ∈U0v U0v = δvi |δvi ∈C 0 , δvi = 0on Γ iv{vi ∈U vU v = vi |vi ∈C 0 , vi = vi on Γ vi}}(test function)(7.13.1a)(trial function)(7.13.1b)To establish the weak form of the momentum equation, Eq.

(7.12.2b), we take theinner product of the momentum equation, Eq. (7.12.2b), with a test function δvi andintegrate over the spatial region, that is:∫Ω δvi ρ∂σ ij∂vi∂vdΩ + ∫ δvi ρc j i dΩ − ∫ δvidΩ − ∫ δvi ρgi dΩ = 0ΩΩΩ∂t χ∂x j∂x j(7.13.2)It is noticed that the stresses in Eq. (7.13.2) are functions of the velocities. Applyingintegration by parts on the stress term yields:∫Ωδvi∂σ ij∂x jdΩ = ∫Ω∂∂ (δvi )(δvi σ ij )dΩ − ∫σ ij dΩΩ∂x j∂x j(7.13.3)By the Gauss divergence theorem, the first term in the RHS of Eq. (7.13.3) can be writtenas a boundary integral, which is:W.K.Liu, Chapter 757∂(δvi σ ij ) dΩ = ∫ t δvi (n j σ ij )dΓΓ∂x j∫Ω(7.13.4)Substituting the specified boundary condition defined in Eq.

(7.12.3b) into Eq. (7.13.4)gives:∫Γtδvi (n j σ ij )dΓ = ∫Γtδvi t j dΓ(7.13.5)The second term in the RHS of Eq. (7.13.3) can be further expressed by using thedefinition of the constitutive equation, Eq. (7.12.5a). That gives:∂(δvi )∂(δvi )σ ij dΩ = ∫− pδ ij + 2µDij dΩΩ∂x j∂x j[∫Ω= −∫Ω∂(δvi )pdΩ+∂xi]∫Ω∂(δvi )2µDij dΩ∂x i(7.13.6)We now may use the decomposition of the velocity gradient in Eq. (7.13.6) into symmetricand antisymmetric parts:∂(δvi )= δLij = δDij + δω ij∂x i(7.13.7a)whereδDij =1  ∂(δvi ) ∂(δv j ) 1  ∂(δvi ) ∂(δv j ) +andδω=−ij2  ∂x j∂xi 2  ∂x j∂xi (7.13.7b)Since ω ij is antisymmetric and Dij is symmetric, it leads to ω ij Dij = 0 .

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

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