© Strojni{ki vestnik 48(2002)12,645-662 © Journal of Mechanical Engineering 48(2002)12,645-662 ISSN 0039-2480 ISSN 0039-2480 UDK 519.61/.64:532.5 UDC 519.61/.64:532.5 Izvirni znanstveni ~lanek (1.01) Original scientific paper (1.01) Metoda robnih elementov za dinamiko viskoelasti~ne Maxwellove teko~ine The Boundary-Element Method for the Dynamics of a Viscoelastic Maxwell Fluid Leopold Škerget - Matej Po`arnik V prispevku je prikazan razvoj numerične sheme na osnovi metode robnih elementov (MRE) za modeliranje ravninskih tokov viskoelastične tekočine. Kot podlaga za razvoj je namenjena shema za modeliranje nestisljivih viskoznih tokov, ki smo jo razširili in ji dodali potrebne člene za zajemanje viskoelastičnosti nenewtonskih tekočin. Posebno pozornost smo namenili integraciji ohranitvenih zakonov in reoloških modelov. Shema je zapisana za hitrostno-vrtinčno formulacijo vodilnih enačb. Kot testni primeri so predstavljeni tokovi nenewtonske Maxwellove tekočine v kanalih različnih geometrijskih oblik. © 2002 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: modeli robnih elementov, dinamika tekočin, tekočine viskoelastične, modeli Maxwellovi) In this paper we show a numerical scheme based on the boundary-element method (BEM) for the numerical modeling of planar, viscoelastic fluid flows. In particular, the singular-boundary-integral approach, which has been established for the viscous, incompressible flow problem, is modified and extended to include the viscoelastic fluid state. Special attention is given to a proper integration of the conservation and constitutive equations. A velocity-vorticity formulation of the governing equations is adopted. As test cases, non-Newtonian-Maxwell fluid flows in the channels of various geometries are studied. © 2002 Journal of Mechanical Engineering. All rights reserved. (Keywords: boundary element methods, fluid dynamics, viscoelastic fluids, Maxwell models) 0 UVOD V prispevku predstavljamo numerični model robnih elementov (MRE) za reševanje dinamike viskoelastične nestisljive tekočine. Sistem Navier-Stokesovih enačb ob upoštevanju Maxwellovih modelov tečenja viskoelastične tekočine smo zapisali za hitrostno-vrtinčno formulacijo. Ta ima v primeru robnoobmočne integralske predstavitve vodilnih enačb določene prednosti, saj smo iz računske sheme izločili tlak, kar pomeni, da se izognemo težavam pri predpisovanju robnih vrednosti tlaka. Posebna pozornost je posvečena spremembi diferencialnih vodilnih enačb v pripadajoče robnoobmočne integralske predstavitve, ki zadoščajo omejitvenemu kriteriju ohranitve mase oziroma pogoju solenoidnosti tokovnega polja. Stabilnost in natančnost numeričnega algoritma dosežemo brez dodatnih stabilizacijskih tehnik z uporabo parabolične difuzijske osnovne rešitve, ki opisuje linearni del prenosnega pojava. Numerični model, predstavljen v prispevku, je mogoče nadgraditi 0 INTRODUCTION This paper deals with a numerical scheme based on the boundary-element method (BEM) developed to numerically simulate viscoelastic, incompressible fluid flows. The method is based on a solution of the Navier– Stokes equations set in the velocity-vorticity formulation for different viscoelastic Maxwell fluids. This formulation has some advantages in the case of a boundary-domain integral representation of governing equations. It is simpler than a primitive variable formulation since the pressure does not appear explicitly in the field-function equations and thus the well-known difficulty connected with the computation of the pressure boundary values in incompressible flows is avoided. Particular attention is given to a proper transformation of the differential governing equations to the corresponding boundary-domain integral representations that satisfy the continuity equation or the solenoidality of the velocity field exactly. The stability and the accuracy of the developed numerical scheme is achieved by employing the parabolic diffusion fundamental solution describing the linear part of transport phenomena without any additional stabilization techniques. The numerical model presented here is easy to upgrade in the sense of effectiveness, gfin^OtJJIMISCSD 02-12 stran 645 |^BSSITIMIGC [kerget L., Po`arnik M.: MRE za dinamiko viskoelasti~ne - BEM for the Dynamics of a Viscoelastic v smislu učinkovitosti, stabilnosti in natančnosti z uporabo difuzivno-konvektivne osnovne rešitve, ki zajema večji del prenosnega pojava, in modela podobmočij [5]. 1 VODILNE ENAČBE 1.1 Ohranitveni zakoni Analitični opis gibanja zvezne snovi temelji na osnovnih ohranitvenih zakonih mase, gibalne količine, energije in kemijskih sestavin, pripadajočih reoloških modelih in enačbah stanja. V obravnavi se bomo omejili na laminarni tok nestisljive viskoelastične izotropne tekočine v območju rešitve W, ograjenim z ograjo G. Opazovane funkcije polja so vektor hitrosti vi(rj,t), tlak p(rj) in temperaturno polje T(rj,t), ki zadoščajo kontinuitetni, gibalni in energijski enačbi: dv dx stability and accuracy using a diffusion-convective fundamental solution where a larger part of the transport phenomena is taken into account. The subdomain technique can also be implemented [5]. 1 GOVERNING EQUATIONS 1.1 Conservation laws The analytical description of the motion of a continuous medium is based on the conservation of mass, momentum, energy and species concentration, with associated rheological models and equations of state. The present development will be focused on the laminar flow of an incompressible, viscoelastic, isotropic fluid in a solution domain W bounded by a boundary G. The field functions of interest are the velocity vector field vi(rj,t), the pressure field p(rj,t), and the temperature field T(rj,t), such that the mass, momentum, and energy equations: = 0 r0 cr p0 0 Dvi = Dt DT j j dx rgi dqj Dt ki so zapisane v obliki kartezičnega tenzorskega zapisa xi in kjer so p0 in cp0 nespremenljiva masna gostota oziroma izobarna specifična toplota, t čas, gi težnostni pospešek, sij napetostni tenzor, qi gostota difuzijskega toplotnega toka in D()/Dt =d()/dt+vkd()/dxk Stokesov oziroma snovski odvod spremenljivke (.). V gibalni enačbi (2) smo upoštevali nestisljivost tekočine v smislu Boussinesqove poenostavitve, pri kateri je vpliv spremenljive masne gostote zajet le v prostorninski sili pgi. Zapisani sistem ohranitvenih enačb pomeni nedokončan sistem parcialnih diferencialnih enačb, ki ga moramo dopolniti z ustreznimi reološkimi modeli oziroma modeli tečenja za posamezno tekočino in znanimi robnimi ter začetnimi pogoji. Ti so v splošnem odvisni od zapisanega sistema enačb in so lahko podani za osnovne fizikalne ali izpeljane funkcije polja. 1.2 Reološki modeli Cauchyjev napetostni tenzor si lahko v primeru toka nestisljive tekočine ločimo na prispevek tlaka in dodatnega deviatoričnega napetostnega tenzorjar: dx (1) (2) (3) j written in the Cartesian frame xi are satisfied where p and cp0 are the constant fluid mass density and the isobaric specific heat capacity, t is the time, gi is the gravitational acceleration vector, sij denotes the components of the total stress tensor, qi stands for the specific heat flux, and D(.)/Dt =3(.)^+vk3(.)/3x represents the Stokes material derivative of the variable (.). The natural convection effect is considered in the momentum, Eq. (2), by using the Boussinesq approximation, where the temperature’s influence on mass density is considered only with the body force pgi. The set of field equations needs to be closed and solved in conjunction with the appropriate rheological equations of the fluid and the boundary, and the initial conditions of the flow problem. The boundary conditions in general depend on the dependent variables applied, i.e. the primitive or velocity-vorticity variables formulation. 1.2 Rheological models For an incompressible fluid the Cauchy total stress si can be decomposed into a pressure contribution plus an extra deviatoric stress-tensor field functions: kjer je Sij Kroneckerjeva funkcija delta. Bistveni del dinamike viskoelastičnih tekočin je izbira primernega reološkega modela, ki podaja soodvisnost dodatnih pdij +tij (4) where dij is the Kronecker delta. The central problem in visco-elastic fluid dynamics is the selection of an appropriate rheological model that relates the extra 2 Sšnn3(aul[M]! ma stran 646 [kerget L., Po`arnik M.: MRE za dinamiko viskoelasti~ne - BEM for the Dynamics of a Viscoelastic napetosti v en. (4) in kinematike toka. Obravnavo bomo omejili na razmeroma preproste diferencialne konstitutivne enačbe, znane kot implicitni Maxwellovi reološki modeli. Za širok spekter snovi, kakor so na primer tekočine s pojemajočim spominom, lahko konstitutivne enačbe zapišemo v obliki odvisnosti med napetostnim in deformacijskim tenzorjem ter njunimi odvodi po času. Za Eulerjev postopek velja, da lahko snovske časovne odvode poljubnega simetričnega tenzorja drugega reda ui , ki ga lahko enačimo z dodatnim napetostnim tenorjem tij ali deformacijskim tenzorjem iij, podamo na različne načine. Z naslednjim izrazom definiramo Stokesov časovni snovski odvod tenzorja: stress in Eq. (4) to the flow kinematics. The differential constitutive equations to be considered here are implicit, rate-type rheological models, generally associated with the name of Maxwell. For the broad class of materials, such as simple fluids with fading memory, the constitutive equation can be expressed using a relation between stress and the strain-rate tensor and their time derivatives. For an Eulerian reference frame the material time derivative or convected derivative of an arbitrary symmetric, second-order tensor uij, which can be equated to the extra-stress tensor tij or the strain-rate tensor e&ij , can be formulated in several ways. Let us first define the Stokes material time derivative with the expression: nato zgornje konvektivni oziroma sodeformabilni: then the upper-convected or codeformation: v Duij u ij = Dt -uik Ljk-ujk Lik (5) (6) while the lower-convected derivative is defined as: (7) medtem ko je spodnje konvektivni odvod podan z izrazom: D Duij uij = Dt + uik Lkj+ujk Lk kjer nadpisa V in D označujeta zgornje oziroma where the superscripts V and D stand for the upper- spodnje konvektivni odvod. Hitrostni gradient L je or lower-convected derivatives, respectively, and the definiran kot: ij velocity-gradient tensor L is defined as: Lij=Lx + Wij (8) kjer je h tenzor vrtilnih hitrosti. Pri implicitnih Maxwellovih reoloskih modelih podamo dodatni napetostni tenzor z vsoto viskoznih in elastičnih učinkov: where W&ij denotes the rotational velocity tensor. Using implicit Maxwell rheological models the extra-stress tensor is given as a sum of the viscosity and elasticity effects: t + t (9) Linearni Maxwellov konstitutivni model podamo z odvisnostjo: The linear Maxwell constitutive model is given as follows: tij=2rl0Lij-\ dt (10) kjer snovska stalnica tekočine 1 pomeni napetostni sprostitveni čas, medtem ko je % dinamična viskoznost. Vpliv elastičnosti, ki je podan z lokalnim časovnim odvodom dodatnega napetostnega tenzorja, je pomemben le med prehodnim pojavom. Z razvojem hitrostnega polja lokalni časovni odvod izgublja pomembnost, tako da v ustaljenem stanju prevladujejo učinki viskoznosti. Najpreprostejši kvazilinearen Maxwellov model, podan v obliki: tij = 2% zajema nelinearnost zaradi lokalnega in konvektivnega odvoda napetostnega tenzorja. where l1 is a material constant for the fluid and is called the stress-relaxation time, whereas h0 stands for the dynamic viscosity. The elasticity is given by the local time derivative of the additional stress tensor, which is significant during transient conditions. As the flow develops, the local time derivative losses its influence to finally arrive at a steady state where the viscosity is dominant. The simplest quasilinear Maxwell model may be given in the following form: Dtij Dt (11) where the nonlinearity of the model is due to the local and convective derivatives of the stress tensor. gfin^OtJJlMlSCSD 02-12 stran 647 |^BSSITIMIGC [kerget L., Po`arnik M.: MRE za dinamiko viskoelasti~ne - BEM for the Dynamics of a Viscoelastic Model predstavljen z en. (11), v praksi ni uporaben, kljub temu pa ga zaradi preprostosti lahko uporabimo pri študiju stabilnosti in natančnosti razvitega numeričnega algoritma. V nadaljevanju sta prikazana oba konvektivna Maxwellova modela. Zgornje konvektivni Maxwellov model podamo z odvisnostjo: The model (11) does not relate to practical problems, but it can be examined due to its simplicity as an appropriate model to study the stability and accuracy of the developed numerical solution algorithm. In addition, both convected Maxwell models are studied. The upper-convected Maxwell model is governed by the following constitutive relation: t 2fl0Lij-\Ty (12) medtem ko za spodnje konvektivni model velja zapis: while for the lower-convected Maxwell model the following equation is valid: D tij = 2h0 e&ij -l1 t ij (13) Zgornje konvektivni Maxwellov model je namenjen testiranju numeričnih modelov reševanja viskoelastičnih tokov, saj je gibanje nekaterih stvarnih tekočin, vsaj v omejenem obsegu, mogoče opisati z en. (12). Preostane še, da zapišemo konstitutivni model difuzijskega toplotnega toka v en. (3). V primeru intenzivnega prehodnega pojava prenosa toplote moramo upoštevati končno hitrost širjenja motnje, tako da velja odvisnost: The upper-convected Maxwell model (12) is used extensively in testing numerical solution models. Some real fluids behave qualitatively like Eq. (12), at least over a limited range of kinematics. Let us write the constitutive model of diffusion heat flux in Eq. (3). In the case of intensive, unsteady heat transfer it is important to take into account the terminal velocity of a moving frontier: dT dx l0 dt (14) kjer sta snovski stalnici k0 in 1 toplotna prevodnost in toplotni sprostitveni čas. Za večino praktično pomembnih prenosnih pojavov zadošča poenostavitev, znana kot Fourierjev zakon difuzije toplote: 1.3 Povzetek vodilnih enačb Z upoštevanjem konstitutivnih modelov za napetostni tenzor (9) do (13) in gostoto toplotnega toka (15) v ohranitvenih zakonih (2) in (3) izpeljemo naslednji sistem nelinearnih enačb: where the material constants k0 and l0 stand for the heat conductivity and the heat relaxation time. For most heat-transfer problems of practical importance the simplification known as the Fourier law of heat diffusion is accurate enough: dT dxi 1.3 Summary of the governing equations (15) Combining the constitutive models for the stress tensor (9) to (13) and the specific heat flux (15) in the conservation equations (2) and (3), the following system of nonlinear equations is developed: dvj dxj = 0 r0 Dvi Dt dp d2vi dti j dx ox ox dx tij = 2h0 e&ij +ti ej dT dt d2T dxjdxj cp0p0 (16) (17) (18) (19) kjer je a0 = k0/ c n toplotna difuzivnost. Sklenjen sistem enačb (16) do (19) je formalno identičen enačbam, ki opisujejo tok newtonske viskozne tekočine z izjemo dodatnega člena, ki je posledica elastičnih učinkov tekočine in je v sistemu označen z nadpisom “e”. Ta člen v enačbi ohranitve gibalne where a0 = k0/ cp0r0 is the thermal diffusivity. The closed system of equations (16) to (19) is formally identical to the equations governing the motion of a Newtonian viscous fluid, except for the additional term describing the elasticity effects (denoted by the superscript “e”). The appearance of the additional terms 2 jgnnatäüllMliBilrSO | | ^SSfiflMlGC | stran 648 [kerget L., Po`arnik M.: MRE za dinamiko viskoelasti~ne - BEM for the Dynamics of a Viscoelastic količine vpelje dodatno nelinearnost v dinamični sistem enačb. Za tok v ravnini sistem enačb od (16) do (19) zagotavlja sedem enačb za sedem neznank, v1, v2, p, t11, t12, t22, in T. Enačbe razrešimo z upoštevanjem primernih robnih in začetnih pogojev. Če predpostavimo, da so vse funkcije polja, npr. , v časovnem koraku t = t znane, je treba določiti vrednosti funkcij polja v naslednjem trenutku tn+1 = tn + Dt 2 HITROSTNO-VRTINČNA FORMULACIJA Z vektorjem vrtinčnosti wi(rj,t) , ki pomeni rotor hitrostnega polja vi(rj): ovk računsko shemo gibanja tekočine razdelimo na kinematični in kinetični del [8]. Kinetiko toka podamo z nelinearno parabolično difuzivno-konvektivno prenosno enačbo vrtinčnosti, ki jo izpeljemo z delovanjem operatorja rotor na gibalno enačbo (17): Dwi dwi dwi d2wi ------ =----- + vj ---- = V0 ------- + wj Dt dt dxj dxjdxj za ij,k,m = 1,2,3. V prenosni enačbi (21) je opazen dodatni prenosni člen wjL , ki ga ne zasledimo v gibalni enačbi (17). Z upoštevanjem odvisnosti: in the momentum equation at the same time increases the nonlinearity of the dynamical system of equations. For the two-dimensional planar geometry the Eqs. (16) to (19) provide seven relations for the seven unknowns, v1, v2, p, t1, t1 t2 and T. The above field equations are to be solved for the appropriate boundary and initial conditions. Assuming that at time level t = tn all the relevant field quantities, such as vn = vn(rj,t), , etc., are known, the issue is to determine the field functions during the time level tn+1 = tn + Dt. 2 THE VELOCITY-VORTICITY FORMULATION Introducing the vorticity vector field function wi(r t) as a curl of the corresponding velocity field vir t): wj = 0 (20) ox j which is a solenoidal vector by definition, the fluid-motion computation procedure is partitioned into its kinetics and kinematics [8] parts. The vorticity transport in the fluid domain is governed by a nonlinear parabolic diffusion-convection equation obtained as a curl of the momentum eq. (17): dvi 1 dp e ijk gh j P0 Sxj p0 dxjdxk for ij,k,m = 1,2,3. In the transport eq. (21) an additional transport term wjL is represented, whereas it cannot be found in eq. (17). Regarding the relation: wjLij=wj(Eij + W& ij) = wjEij (22) preprosto pokažemo, da del člena w J&ij, i&ij za i = j, opravlja prenos vrtinčnosti z raztezanjem vrtinčne črte, preostali del, i&ij za i*j, pa prenaša vrtinčnost z obračanjem vrtinčne črte. Prenosni pojav raztezanja in obračanja ne obstaja pri ravninskih tokovih, ko se vektorska enačba (21) poenostavi v skalarno: Dw =dw dw = d w Dtdt +vjdxV0dxdx it is easy to understand the behaviour of w ji&ij, i&ij for i = j, representing vorticity transport due to the stretching of the vorticity line, while the other part, e&ij for i *j, stands for the twisting of the vorticity line. Twisting and stretching transport phenomena do not exist during the planar flows when vector eq. (21) is reduced to the scalar one: (23) za ij,k = 1,2. Enačba (21) podaja časovno spremembo vrtinčnosti delca tekočine zaradi viskozne difuzije, konvekcije, vzgona, učinkov deformacije in elastičnosti tekočine. V primerjavi z enačbo gibanja newtonske viskozne tekočine je enačba (21) močneje nelinearna, prav zaradi močno nelinearnega elastičnega izvirnega člena. Z delovanjem operatorja rotor na vektor vrtinčnosti: 1 S/> 1 d2tej k +-----eg— +-----e------ P0 ij j 8xi p0 ij dxidxk for ij,k = 1,2. The vorticity equation (21) expresses time-dependent vorticity transport by viscous diffusion, convection, bouyancy forces, while the elasticity and deformation of the fluid acts as a highly nonlinear production term, making the nonlinearity of the equation even more severe, when compared to Newtonian viscous fluid flow. By applying the curl operator to the vorticity definition: Vxw = Vx(Vxv)=V(V-v)-Dv (24) in ob upoštevanju kontinuitetne enačbe (16), div v = 0, izpeljemo naslednjo vektorsko eliptično Poissonovo enačbo: and by using the continuity eq. (16), div v = 0 , the following vector elliptic Poisson equation is obtained: stran 649 [kerget L., Po`arnik M.: MRE za dinamiko viskoelasti~ne - BEM for the Dynamics of a Viscoelastic d2vi ------- + dxjdxj Enačba podaja kinematiko toka nestisljive tekočine oziroma združljivosti in omejitveni pogoj med solenoidnima vektorskima poljema hitrosti in vrtinčnosti v dani točki prostora in časa. Z namenom, da bi povečali konvergenco in stabilnosti vezane hitrostno-vrtinčne iterativne sheme, en. (25) zapišemo v njeni nepravi prehodni obliki, kar se izraža v naslednji parabolični difuzijski enačbi za vektor hitrosti ijk dx (25) The equation represents the kinematics of an incompressible fluid motion expressing the compatibility and restriction conditions between the solenoidal velocity and the vorticity vector field functions at a given point in space and time. To accelerate convergence of the coupled velocity-vorticity iterative scheme the false transient approach is used. Thus, in a solution scheme Eq. (25) is rewritten as parabolic diffusion equation for the velocity vector: d v 1 dv dak dxjdxj a dt dx (26) kjer je a sprostitveni parameter. Očitno je, da je vodilna with a as a relaxation parameter. It is obvious that the hitrostna enačba (26) natančno izpolnjena v governing velocity equation (26) is exactly satisfied ustaljenem stanju (t-*»), ko umetni časovni odvod only during the steady state (t-*»), when the time hitrosti odpade. derivative vanishes. 3 TLAČNA ENAČBA Za nestisljive tekočine velja ugotovitev, da je tlak le funkcija hitrostnega polja in nasprotno, medtem ko je gostota nespremenljiva. Tlak torej ni termodinamična veličina stanja. Zapišimo gibalno en. (17) z upoštevanjem vektorske enakosti (24) za tlačni gradient grad p: dxi fi P0 Dt ki se v primeru ravninskega toka poenostavi v odvisnost: 3 PRESSURE EQUATION In the case of incompressible fluid motion, pressure is only a function of the velocity field, and vice versa, while the fluid mass density is assumed to be invariant. Therefore, pressure is not a thermodynamic variable. Let us write momentum, Eq. (17), for the pressure gradient, grad p, by taking into account the vector relation (24): dmk dri j dx dx (27) which simplifies in the case of planar flows in the following way: dp dx fi = -r0 Dvi Dt da) dri j ?l0eij — + Pgi+— ox ox (28) V vektorski funkciji fi smo združili vztrajnostne, difuzijske, težnostne in elastične učinke. Pri izpeljavi tlačne enačbe v odvisnosti od hitrosti poiščemo divergenco členov en. (27), kar se kaže v naslednji eliptični Poissonovi enačbi za tlak: Inertia, diffusion, gravitational and elasticity effects are incorporated in vector function fi. To derive the pressure equation dependent on velocity the divergence of Eq. (27) should be calculated, resulting in the elliptic Poisson pressure equation dp dfi dxdx dx (29) 4 METODA ROBNIH ELEMENTOV 4 THE BOUNDARY-ELEMENT METHOD Uporaba Greenovih osnovnih rešitev kot utežnih funkcij v numeričnem modelu pomeni eno izmed osnovnih lastnosti in prednosti metode robnih elementov (MRE). Ker osnovne rešitve opisujejo le linearni del prenosnega pojava, je izbira primernega linearnega operatorja L[.] ključnega pomena za zapis integralskih enačb, ki ustrezajo prvotnemu sistemu ohranitvenih diferencialnih enačb. Znane diferencialne modele prenosnih pojavov v toku tekočine lahko zapišemo v obliki naslednjega splošnega stavka: The unique advantage of the boundary-element method (BEM) originates from the application of Green fundamental solutions as particular weighting functions. Since they only consider the linear transport phenomenon, the appropriate selection of a linear differential operator L[.] is of vital importance in establishing stable and accurate singular integral representations corresponding to original differential conservation equations. All the various conservation models can be written in the following general form L[u]+b =0 (30) 2 jgnnatäüllMliBilrSO | | ^SSfiflMlGC | stran 650 [kerget L., Po`arnik M.: MRE za dinamiko viskoelasti~ne - BEM for the Dynamics of a Viscoelastic kjer je L[.] eliptični ali parabolični linearni diferencialni operator, u(rj) je poljubna funkcija polja in b(r,t) pomeni nehomogeni del nelinearnih vplivov oziroma prostorninskih virov. 4.1 Integralska predstavitev kinematike where the operator L[.] can be either elliptic or parabolic, u(rj,t) is an arbitrary field function, and the nonhomogenous term b(rj,t) is applied for nonlinear transport effects or pseudo body forces. 4.1 Integral representation of flow kinematics Hitrostno en. (26) lahko prepoznamo kot The velocity Eq. (26) can be recognized as a parabolično difuzijsko enačbo, zato uporabimo parabolic diffusion equation. Employing the linear linearni parabolični difuzijski operator: parabolic diffusion operator as follows: Lfl = aP) -?L) tako, da lahko zapišemo: dxjdxj dt the following can be stated: Lv]+b = aA2vi -^+b = 0 i j i dxjdxj dt (31) (32) Pri izpeljavi singularne robne integralske The singular boundary integral representation enačbe izhajamo iz Greenovih teoremov za skalarna for the velocity vector can be formulated by using polja ali iz integrala utežnega ostanka, kar se kaže v the Green theorems for scalar functions or weighting naslednji obliki vektorske integralske formulacije: residuals technique rendering the following vector integral formulation: (x)vi( x,tF) + a\ vi------dtdG = a\ iu*dtdG+\ bu*dtdW+ \ viF1u*F1dW dn (33) kjer je u* parabolična difuzijska osnovna rešitev, npr. with u* the parabolic diffusion fundamental solution rešitev enačbe: as the solution of the equation: 8 u* 8u* a +-----+ d (x,s)d (tF,t) = 0 dxjdxj dt in podana z izrazom: and given by the expression: u* (4p at ) 1 | r d/2exp 4at (34) (35) kjer sta (x,tF) in (s,t) izvirna točka in referenčna točka where (x,tF) and (s,t) are the source point and the polja, d pomeni izmero problema in t = t - t. Če reference field point, d is the dimension of the problem izenačimo prostorninske sile z enakostjo: and t = tF - t. The pseudo-body force vector includes the vortical fluid flow part: b = ae ijk dx (36) velja končna oblika integralske predstavitve rendering the following final integral formulation of kinematike toka, npr. v vektorski obliki: flow kinematics, in vector form: c ( x ) v ( x,tF ) + a\ v-----dtdG = a\u*dtdG-a\ wxnu*dtdG dn GtF-1 Qn +a\ wx'Vu*dtdW+ vF-1uF-1dW oziroma v tenzorskem simbolnem zapisu W tF-1 ^ or in tensor symbolic notation: c(x)vi ( x,tF) + ajX viJudtdG = al\tFi iu*dtdG du* -aeijk\ \F w jnku*dtdG+ aeijk \ wj-----dtdW+ \ vi, F -1uF -1 dW (37) (38) dx Kinematika ravninskega toka je podana z The kinematics of plane motion is given by dvema skalarnima enačbama: two scalar equations as follows: c(x)vi (x,tF) + a\G f vi udtdG = a\G f viu*dtdG dn G tf-1 dn +aeij w n ju* dtdG - aeij\ w------dtdW+ \ vi F1u*F1dW (39) stran 651 glTMDDC [kerget L., Po`arnik M.: MRE za dinamiko viskoelasti~ne - BEM for the Dynamics of a Viscoelastic Eno osnovnih izhodišč numeričnega modela toka nestisljive tekočine je divergence prosta oziroma solenoidna rešitev za hitrostno in vrtinčno polje. Za en. (36) oziroma (38) lahko preprosto pokažemo, da sta enačbi izpolnjeni tudi v primeru, ko nobena od divergenc ni enaka nič [9]. Povzamemo lahko, da en. (38) v splošnem ne predstavlja kinematike toka nestisljive tekočine. S ponovno zahtevo po združljivosti in omejitvi hitrostnega ter vrtinčnega polja, r> = rotv in divv = 0, robna integrala na desni strani en. (38) preoblikujemo, kar se kaže v končni obliki singularne robne integralske predstavitve kinematike prostorskega toka: One of the most important issues in the numerical modeling of incompressible flow phenomena is to obtain a divergence-free final solution both for the velocity field and for the vorticity conservation. In the case of equations (36) or (38) it can be easily shown that solutions where none of the divergences are zero are permitted [9]. Thus, it is possible to conclude that Eq. (38) does not generally represent the kinematics of incompressible fluid flow. By using additional compatibility and restriction conditions for the velocity and vorticiy fields, w = rot v and div v=0, boundary integrals on the right-hand side of Eq. (38) are rewritten, resulting in the final singular boundary integral statement of the kinematics of spatial flow: (x)v(x,tF) + a\ \F (Vu* -n ) vdtdG = a\ \tF (Vu* xn )x vdtdG + a\ \tF \ w x Vu* J dtdW + j vF-u*-dW (40) ali tudi v zgoščenem simbolnem zapisu za krožno kombinacijo indeksov ijkij = 12312: or in the compact symbol notation for the cyclic combination of: ijkij = 12312: c(x)vi ( x,tF) + a\XviudtdG = al\t; du* du* vk -------ni---------nk 1 { Öxk dxi dtdG du du vj ------nj--------ni 1 <9x dx dtdG + a\ \F wj-----dtdW-a\ \F co-----dtdW + \ JW lt JW lt k JW (41) v,- p^u^dW -a\ \tF JGJtF -1 j[dxi j dxj i\ WtF -1 j dxk WtF-1 k dxj Kinematika ravninskega toka je podana z: The kinematics of planar fluid motion is given by: c(x)vi ( x,tF) + a\G\tF vi^dtdG = a\G\tF vudtdG-ae^JtF wudtdW + jW vi,F -1uF -1dW (42) dn Enačba (41) je popolnoma ustrezna kontinuitetni oziroma omejitveni enačbi in definiciji vrtinčnosti ter podaja kinematiko nestisljive tekočine v integralski obliki. Hitrostni robni pogoji so vključeni v robnih integralih, medtem ko je z območnim integralom zajet vpliv vrtinčnega polja na razvoj hitrostnega polja. Zadnji območni integral upošteva vpliv začetnih hitrostnih pogojev nepravega prehodnega pojava. V en. (42) se pojavljata normalni in tangentni odvod osnovne rešitve du*/dn in du*/dt, kar je pomembna razlika proti en. (39), kjer sta osnovna rešitev in njen normalni odvod, u* in du*/dn. Pri izračunu robnih vrednosti funkcij polja moramo uporabiti normalno oziroma tangentno obliko vektorske enačbe (40), [7]. Robne vrednosti vrtinčnosti so v integralski obliki zajete v območnem integralu. Izračun robnih vrednosti vrtinčnosti zaradi zapisa nesingularnega implicitnega sistema terja tangentno obliko vektorske enačbe (40). dxj Equation (41) is equivalent to the continuity equation, also recognized as the restriction equation and the vorticity definition expressing the kinematics of general incompressible fluid flow in the integral form. Velocity boundary conditions are incorporated into boundary integrals, while in domain integrals the influence of the vorticity field on the developing velocity field is given. The last domain integral takes into account the influence of the initial velocity conditions of false transient phenomena. In Eq. (42) the normal and tangential derivatives of the fundamental solution,5u*/cn and5u*^, are employed an important difference in comparison to Eq. (39), where the fundamental solution and the normal derivative of the fundamental solution, u* and 3u*/3n, are used. To compute the boundary values of the field functions, the normal or tangential form of vector Eq (40), [7], is required. The boundary vorticity values are expressed in integral form within the domain integral. When the unknows are the boundary vorticities one has to use the tangential component of vector Eq. (40) because of the nonsingular implicit system of equations c(x)n(x)xv(x,tF) + n(x)xaj jtF (Vu* -n)vdtdG = n ( x )x.a\ V (Vu*xn)xvdtdG + n (^)xa\ \F lwxVu*\dtdW + nx )xa\ vF-1uF -1dW 4.2 Integralska predstavitev kinetike 4.2 Integral representation of flow kinetics (43) Za zapis kinetike toka viskoelastične tekočine Considering the kinetics of viscoelastic fluid v integralski obliki moramo upoštevati parabolično flow in an integral representation one has to take into 2 jgnnatäüllMliBilrSO | | ^SSfiflMffiC | stran 652 [kerget L., Po`arnik M.: MRE za dinamiko viskoelasti~ne - BEM for the Dynamics of a Viscoelastic difuzivno-konvektivno naravo prenosne enačbe vrtinčnosti. Z uporabo linearnega paraboličnega difuzivnega diferencialnega operatorja: L[-] prenosno enačbo vrtinčnosti (33) zapišemo v obliki nehomogene parabolične difuzivne enačbe: account the parabolic diffusion-convection character of the vorticity transport equation. Since only the linear parabolic diffusion differential operator is employed, i.e.: dxjdxj dt (44) the vorticity equation (33) can be formulated as a nonhomogenous parabolic diffusion equation as follows: d2w dw L[w] + b = v—--------— + b = 0 dxdx dt (45) z naslednjim pripadajočim integralskim stavkom, zapisanim za časovni korak Dt = tF - tF-1: with the following corresponding integral representation written in a time-increment form for a time step Dt = tF - tF-: (č)w (č,tF) + v0\ w-----dtdG = V0\ —ukdtdG+\ bu*dtdW+ \ wF- 1u*F - 1dW -1 dn (46) kjer je u* difuzivna osnovna rešitev, podana z en. (35). Območni integral nehomogenega nelinearnega prispevka b: dvjw 1 b =-^ + — eijgj dxj r0 vsebuje konvekcijo, vzgonske in elastične učinke, tako da velja naslednji zapis kinetike vrtinčnosti v integralski obliki: where u* is the parabolic diffusion fundamental solution given by Eq. (35). The domain integral of the non-homogenous nonlinear contribution b, represented as: sv dr 1 dxi r0 dxidxk (47) includes the convection, the bouyancy force effects, and the viscoelastic effects. Thus the final integral statement reads as: ( ) ( ) r rtF du* r rtF dw r rtF r rtF du* nffltu + yj w----dtdG = vA —u*dtdG-\ wvnu*dtdG +\ wvj----dtdW JG Jt JG Jt JG h JW h F-1 9n +eij — i itF nigjru*dtdG-eij — [ \ gjr-----dtdW + eij — [ [tF 54 cx *dtdG (48) -e— —--------dtdW+\ wF1uF1dW JW JtF-1 JW r0W tF -1 dxk dxi Iz en. (48) je razvidna popolna podobnost med Eq. (48) shows the analogy between the prenosom vrtinčnosti v viskoelastični tekočini in vorticity transport in the viscoelastic fluid and the prenosom vrtinčnosti v toku newtonske viskozne vorticity transport in viscous, Newtonian motion, tekočine z izjemo dodatnega elastičnega prispevka, with the only difference in the extra viscoelastic ki deluje kot močno nelinearni izvirni člen. contribution acting as a highly nonlinear source term. Singularno robno integralsko predstavitev By applying a similar procedure to the heat- toplotne prenosne enačbe izpeljemo enako kakor izpeljavo transport equation, one derives the following integral vrtinčne enačbe, tako da velja integralski stavek: statement: c(^)T(^,tF) + a0jJtF T—dtdG = a0jJtF —u*dtdG-jGjF Tvnu*dtdG dn G tf-1 dn + { [tF Tvj—dtdW+\ TF1u*F1dW JW hF1 JW (49) 4.3 Integralska predstavitev tlačne enačbe 4.3 Integral representation of the pressure equation Tlačna en. (29) je eliptična Poissonova enačba, The pressure Eq. (29) is recognized as an tako da uporabimo linearni eliptični Laplacev elliptic Poisson equation, thus employing the linear diferencialni operator: elliptic Laplace differential operator: L[] kar se izraža v zapisu: g2 () dxidxi the following can be stated: LU + b = Üp + b = 0 [] dxdx (50) (51) ^vmskmsmm 02-12 stran 653 |^BSSITIMIGC [kerget L., Po`arnik M.: MRE za dinamiko viskoelasti~ne - BEM for the Dynamics of a Viscoelastic medtem ko je pripadajoča singularna robna integralska tlačna predstavitev podana s stavkom: while the corresponding singular integral pressure representation is given by: du* dp () p ($) +\Gp—dG = lfn u*dG+lbu*dW (52) kjer je u* Laplaceova osnovna rešitev, npr. rešitev enačbe: where u* is the Laplace fundamental solution, being the solution of the equation: dxdx + S(š,s (53) Za ravninsko geometrijo velja rešitev: Z izenačitvijo navidezno prostorninskih sil: 2n For planar geometry the following solution is suitable: (54) By equating pseudo-body forces with the expression: Ofi