- Email: [email protected]
RBF meshless modeling of non-Newtonian Hele–Shaw flow
ARTICLE IN PRESS
Engineering Analysis with Boundary Elements 31 (2007) 863–874 www.elsevier.com/locate/enganabound
RBF meshless modeling of non-Newtonian Hele–Shaw flow Francisco Bernal, Manuel Kindelan Universidad Carlos III de Madrid, Avda. Universidad 30, 28911 Leganes, Spain Received 31 March 2006; accepted 31 January 2007 Available online 20 April 2007
Abstract In this paper we consider the problem of injecting a non-Newtonian fluid into a thin cavity. Using the Hele–Shaw approximation the problem reduces to a moving boundary problem in which the pressure is described by a 2D, nonlinear, elliptic equation. Mesh-free methods are very well suited for the numerical solution of moving boundary problems since no remeshing is needed at each time step to correctly represent the boundary. Among these methods, we have chosen the asymmetric RBF collocation method (Kansa’s method) to compute the pressure distribution. Kansa’s method is truly mesh-free, and is one of the most frequently used methods due to its accuracy and ease of implementation. Once the pressure is known, the velocity at each point in the moving front is computed and, therefore, the location of the front can be updated. To efficiently model the front motion we use the level set method, which is very accurate, and can handle both front collisions and front break-ups without difficulty. r 2007 Elsevier Ltd. All rights reserved. Keywords: Radial basis function; Generalized multiquadries; Level set method
1. Introduction Recently, very intensive efforts have been devoted to develop mesh-free or element-free methods that eliminate the need of element connectivity in the solution of partial differential equations (PDEs). The motivation is to cut down modeling costs in industrial applications by avoiding the labor intensive step of mesh generation. These methods are particularly attractive in problems with moving interfaces since no remeshing is necessary. In this paper we use the radial basis function (RBF) method which was first used by Hardy [1] for interpolation of scattered data and later by Kansa [2,3] for the solution of PDEs. Kansa’s method looks for an approximate solution in the function space spanned by a finite set of RBFs, by using collocation of the PDE in a set of (scattered) centers. In particular, we use multiquadric (MQ) RBFs which have been shown to have exponential convergence [4,5]. However, the price that one has to pay for this increased accuracy is ill-conditioning of the resulting system of Corresponding author.
E-mail address: [email protected] (M. Kindelan). 0955-7997/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2007.01.008
algebraic equations which have to be solved to determine the coordinates of the approximate solution in the radial functions basis. In fact, there is a trade-off between accuracy and ill-conditioning which Schaback explains by a kind of uncertainty relation [6]. During the last years many procedures have been proposed to improve the conditioning problem; approximate cardinal basis function (ACBF) preconditioning [7–10], nonoverlapping [11] and overlapping [12] domain decomposition methods, differential quadrature methods [13], affine space decomposition [14]. The fact is that, although the RBF method has many attractive properties, it has not been yet extensively applied to solve practical problems. Even more, most of the engineering problems [15–20] where the method has been used are linear ones. Recently, Fasshauer [21] showed that Kansa’s method can be easily adapted for nonlinear elliptic PDEs. In this paper we address the problem of injection molding, which is one of the most important industrial processes for the manufacturing of thin, plastic products. The constant demand for higher quality leads to a growing interest in the analysis and prediction of the product’s physical properties through the use of numerical simulations.
ARTICLE IN PRESS 864
F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
The mathematical model of injection molding is a free boundary problem described by conservation equations for mass, momentum and energy. The fluid is assumed to be non-Newtonian and a power-law equation is used to model its rheology. The conservation equations for this problem, can be dramatically simplified assuming that the mould is thin so that a Hele–Shaw approximation [22,23] can be used. In this case the momentum equation is just a 2D, nonlinear, elliptic equation whose solution yields the pressure distribution in the filled region of the mould. From this pressure field, the velocity distribution can be computed and the location of the advancing front can be updated. The front motion is modeled using a level set method [24], which has proved to be very accurate, and which is able to handle both front collisions and front breakups without difficulty. It is based in considering the moving front as the zero level set of a higher dimensional function which is evolved using the advection equation. For simplicity, the advection equation is numerically solved using the upwind method in a regular equispaced grid. 2. Model equations To model the flow of an injected fluid in a mould, we start from the conservation equations for mass, momentum and energy, qr þ r ðrvÞ ¼ 0, qt qu qu qu qu þu þv þw r qt qx qy qz qp qsx qtxy qtxz þ ¼ þ þ , qx qx qy qz qv qv qv qv þu þv þw r qt qx qy qz qtxy qsy qtyz qp þ ¼ þ þ , qy qx qy qz qw qw qw qw þu þv þw r qt qx qy qz qp qtxz qtyz qsz þ ¼ þ þ , qz qx qy qz qT qT qT qT þu þv þw rC p qt qx qy qz 2 q T q2 T q2 T þ 2 þ 2 þ ZF, ¼k qx2 qy qz
ð2Þ ð3Þ
where the inertia terms in the momentum equations have been neglected in comparison to the viscous terms (creeping flow), and where partial derivatives in the horizontal directions (x; y) have been neglected in comparison to derivatives in the vertical direction (small thickness). Integrating the momentum equations along the vertical direction from the mid-plane, and taking into account that by symmetry ðqu=qzÞðz ¼ 0Þ ¼ ðqv=qzÞðz ¼ 0Þ ¼ 0, results in qu z qv z ¼ Lx ; ¼ Ly , qz Z qz Z qp qp Lx ; Ly . qx qy A new integration in the vertical direction yields, Z b Z b Z b qu Lx z~ z~ d~z ¼ d~z ) u ¼ Lx d~z; q~ z Z Z z z z Z b Z b Z b Ly z~ qv z~ d~z ¼ d~z; d~z ) v ¼ Ly q~ z Z Z z z z
ð4Þ ð5Þ
where bðx; yÞ is the half-height and where the no slip boundary condition, uðz ¼ bÞ ¼ vðz ¼ bÞ ¼ 0, has been used. Integrating these equations along the vertical direction results in the following averaged horizontal velocities: Z Z Ly 1 b Lx 1 b u¯ ¼ S; v¯ ¼ S, (6) u dz ¼ v dz ¼ b 0 b 0 b b where Z S¼ 0
b 2
z dz. Z
Thus, the velocity vector is in the direction of rp. Analogously, the continuity equation (1) is integrated along the half-height, Z b qu qv qðb¯uÞ qðb¯vÞ þ þ . (7) dz ¼ qx qy qx qy 0
where F is the dissipation function. Under the assumptions of incompressibility, creeping flow (ðr1 U 1 L=Z1 Þ ðh=LÞ2 51) and small thickness (h=L51), the above equations simplify to, qu qv þ ¼ 0, qx qy
qp q qu qp q qv qp ¼ Z ¼ Z ¼ 0, ; ; qx qz qz qy qz qz qz qT qT qT q2 T rC p þu þv ¼ k 2 þ ZF, qt qx qy qz
ð1Þ
Introducing the values of the average horizontal velocities (6) into the integral form of the continuity equation (7) produces, q qp q qp S S þ ¼ 0. (8) qx qx qy qy This is an elliptic, nonlinear equation whose solution yields the pressure field and from (4)–(5) the velocity field. Notice that, in general, the viscosity is a function of pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi shear stress (_g ¼ u2z þ v2z ) and temperature (Z ¼ ZðT; g_ Þ).
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
From (4)–(5) g_ ¼
Lz Z
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where L jrpj ¼ L2x þ L2y .
(9)
For instance, if Z is given by the power-law model, namely Z ¼ Z0 gðTÞ_gn1
(11)
In the general case, Z depends on temperature and, therefore, the momentum equation (8) is coupled to the energy equation (3). However, in what follows we will neglect the temperature dependence thus decoupling energy and momentum equations. Taking gðTÞ ¼ 1 the flowability is simply given by nbð2 nþ1Þ=n 1=n Z0 ð2 n
.
(12)
þ 1Þ
We will also assume that the mould width is constant (bðx; yÞ ¼ const.), so that Eq. (8) simplifies to, q qp q qp jrpjm jrpjm þ ¼ 0, (13) qx qx qy qy where m ¼ ð1 nÞ=n. The case m ¼ 0 corresponds to a Newtonian fluid and the case ma0 to a non-Newtonian fluid. From (6) the averaged horizontal velocities are given by u¯ ¼ jrpjm
qp ; qx
v¯ ¼ jrpjm
qp . qy
(14)
Eq. (13) has to be solved with the following boundary conditions:
Many different RBF basis functions have been used in the past. Among them, we have chosen the MQ basis functions, fk ðxÞ ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kx xk k2 þ c2k
inlet pressure (or flow rate) at the mould gate, p ¼ 0 in advancing front, normal velocity zero at the walls of the mould, qp=qn ¼ 0.
The solution of Eq. (13) is used in (14) to compute the average velocity field, which is in turn used to advance the front. This model is the Hele–Shaw approximation for injection molding which was introduced by Hieber and Shen [22,23] and later used by many authors (see for instance [25–27]).
first proposed by Hardy [1] and which have been used extensively both for interpolation and for the solution of PDEs. In a thorough study of interpolation methods for 2D scattered data, Franke [28] found that MQ generally performs best. To specify the MQ basis (16), one must specify both the centers xk and the shape parameters ck . The problem of choosing the optimal shape parameters has been investigated by several authors, although there is no yet a widely accepted theory for choosing them. Recently, Wertz et al. [29] reviewed those efforts and carried out many numerical experiments to study the role of generalized MQ shape parameters in the solution of elliptic PDEs. They conclude that the shape parameters for RBF centers located in the boundary should be much greater than the shape parameters in interior RBF centers. Having specified the MQ basis functions, the problem reduces to determining the unknown coefficients ak ðk ¼ 1; . . . ; NÞ. They are computed by collocation in N ¯ ¼ fx¯ 1 ; x¯ 2 ; . . . ; x¯ N g in the domain, O, collocation nodes, X and in its boundary, qO, where either the PDE (13) or the appropriate boundary conditions are imposed. Usually, the collocation nodes, x¯ i and the RBF centers xi coincide, although it is possible to use different collocation nodes and RBF centers (see for instance [30]). For convenience, Eq. (13) can be rewritten as jrpjm Dp þ
qjrpjm qp qjrpjm qp þ ¼ 0, qx qx qy qy
3. Pressure field To solve Eq. (13) we use an RBF based method. The idea is to find an approximate solution in the function space S ¼ spanffðk xk kÞ : xk 2 Xg, where f is a globally supported RBF and X a set on N RBF centers. Thus, N X k¼1
ak fk ðxÞ;
fk ðxÞ fðkx xk kÞ.
(15)
(17)
which for Newtonian fluids (m ¼ 0) reduces to the Laplace equation. Substituting (15) into (17) and using collocation at node xi , results in, jrpi jm Dpi þ
qjrpi jm qpi qjrpi jm qpi þ ¼ 0, qx qx qy qy
where pi ¼
N X
ak fk ðx¯ i Þ,
k¼1
pðxÞ ¼
(16)
(10)
the flowability S becomes [23] Z 1 L 1=n b z 1=n z dz. S¼ L Z0 g 0
S ¼ Lð1nÞ=n
865
N N qpi X qfk qpi X qfk ¼ ðx¯ i Þ; ¼ ðx¯ i Þ, ak ak qx qx qy qy k¼1 k¼1 2 N X q fk q2 f k ak ð x Þ þ ð x Þ , Dpi ¼ ¯ ¯ i i qx2 qy2 k¼1 2 !2 !2 3m=2 N N X X qf qf ik ik 5 , jrpi jm ¼ 4 ak þ ak qx qy k¼1 k¼1
(18)
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
866
2 !2 !2 3m=21 N N X X qjrpi jm qf qf ik ik 5 ¼ m4 ak þ ak qx qx qy k¼1 k¼1 " ! ! N N X X qfik q2 fik ak ak qx qx2 k¼1 k¼1 ! !# N N X X qfik q2 fik þ ak ak , qy qyqx k¼1 k¼1 2 !2 !2 3m=21 N N X X qjrpi jm qf qf ik ik 5 ¼ m4 ak þ ak qy qx qy k¼1 k¼1 " ! ! N N X X qfik q2 fik ak ak qx qxqy k¼1 k¼1 !# ! N N X X qfik q2 fik þ ak ak , qy qy2 k¼1 k¼1 fi;k ¼ fk ðx¯ i Þ;
qfik qfk ðxÞ ðx¯ i Þ; ¼ qx qx
q2 fik q2 fk ðxÞ ¼ ðx¯ i Þ; qx2 qx2
qfik qfk ðxÞ ðx¯ i Þ, ¼ qy qy
q2 fik q2 fk ðxÞ ¼ ðx¯ i Þ, qy2 qy2
q2 fik q2 fk ðxÞ ðx¯ i Þ, ¼ qxqy qxqy (see the Appendix). Thus, Eq. (18) is a nonlinear algebraic equation on the N unknowns ak , which is imposed in N p PDE collocation nodes x¯ i ; i ¼ 1; . . . ; N p . In addition, there are N b algebraic equations resulting from collocation of the boundary conditions on N b boundary nodes x¯ i ; i ¼ 1; . . . ; N b . If N p þ N b ¼ N, the problem reduces to the solution of N nonlinear equations on the N unknowns ak ; k ¼ 1; . . . ; N. 3.1. Numerical method To solve Eq. (18), we use an iterative procedure based on linearizing equation (18) according to, ðnþ1Þ m jrpðnÞ þ i j Dpi
ðnþ1Þ ðnþ1Þ m m qjrpðnÞ qjrpðnÞ i j qpi i j qpi þ ¼ 0, qx qx qy qy (19)
where the superscriptðnÞ refers to the iteration number. Eq. (19) represents a system of N p linear equations on the N unknowns ak ; k ¼ 1; . . . ; N. A similar procedure is used to linearize the equations resulting from Neumann boundary conditions and, therefore, they provide an additional set of N b linear equations. Thus, at each iteration a system of N linear equations has to be solved to compute the N RBF coefficients ak . We consider the rectangular mould shown in Fig. 1, where the injected fluid is shown in dark. At the inlet gate we specify either the pressure (p ¼ 1, Dirichlet) or the normal derivative (qp=qx ¼ f ðyÞ, Neumann). In the moving
p=0
p=1 ∂ p/∂ x=f(y)
∂ p/∂ x=0
Fig. 1. Problem geometry.
front we take p ¼ 0 (Dirichlet) and at the walls we assume that the normal velocity is zero (qp=qn ¼ 0). For the fluid we define m ¼ 0:6 which corresponds to a non-Newtonian fluid with power-law exponent n ¼ 0:625. In a previous paper [31] dealing with the Hele–Shaw problem in the case of Newtonian fluids, we analyzed this problem using PDE collocation at interior nodes and boundary condition collocation at boundary nodes. We found that large errors occur near the mould entrance, where there is a discontinuity in the boundary condition. These inaccuracies near the boundaries is a well-known feature [30] in all RBF approximations. An alternative to improve the accuracy of the method was proposed by Fedoseyev et al. [32]. It is based on the observation that the residual error is typical largest near the boundaries. Thus, they proposed to enforce collocation of the PDE in boundary nodes so that both the boundary condition and the PDE are imposed in those nodes and the residuals dramatically decrease. However, since the number of equations increases, it is necessary to introduce additional RBF centers to match the number of unknown coefficients ak . In [31] we showed that using this method, or the Not-a-Knot method proposed by Fornberg [30], results in great improvements in accuracy. Fig. 2 shows the initial distribution of collocation nodes and RBF centers. There are N o ¼ 50 interior collocation nodes (red circles), N b ¼ 28 boundary collocation nodes (blue and black circles), and N ¼ 106 RBF centers (crosses). The interior and boundary nodes are used both as collocation nodes and as RBF centers. However, for each collocation node in the boundary, a new RBF center is introduced which is located in the direction normal to the boundary at a distance equal to the average distance between centers. PDE collocation is enforced in N p ¼ N o þ N b ¼ 78 nodes, while boundary condition collocation is enforced in N b ¼ 28 nodes. Thus, there is a total of N ¼ 106 nonlinear equations for the N ¼ 106 unknowns ak ; k ¼ 1; . . . ; N. The N p PDE collocation nodes are taken from a finite element mesh.
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
To start the iterations we use as a first approximation the solution of the Laplace equation (m ¼ 0). Then we proceed with the iterations checking the convergence of the procedure by computing the residual at the collocation nodes, Rðnþ1Þ ðx¯ i Þ ¼ jrpðnþ1Þ jm Dpðnþ1Þ þ þ
resulting from RBF collocation is solved faster than the lower shape parameter of the RBF basis functions. But does this mean that the error of the solution decreases with decreasing shape parameter? To check this, we use Eq. (20) to compute the residual, not at the collocation nodes, x¯ i , but at the nodes of a fine, equidistant mesh. If the computed RBF solution is correct, the residual should be zero at Rany point. To compare the accuracy of the solution we use O jRðxÞj dx which is shown in the right of Fig. 3 as a function of the shape parameter. Notice that the residual decreases with increasing shape parameter until a minimum is reached at approximately c ¼ 0:04. If the number of RBF centers increases, the value of the optimal shape parameter decreases and the residual decreases. This is shown by the dashed line which corresponds to N ¼ 337 RBF centers. A similar behavior was observed in the linear case (m ¼ 0, see [31]). The left part of Fig. 4 shows the RBF solution (15) on a fine, regular mesh. The right part shows the corresponding residuals and the location of the RBF centers. Notice that the residuals are different from zero only in the entrance region, where the discontinuity in the second derivatives cannot be correctly reproduced in the functional space spanned by the smooth RBF basis functions. Notice also, that the residuals are zero at the RBF centers, since they are also used as collocation nodes, but they exhibit large oscillations between them. Fig. 5 shows the first derivatives of the RBF solution (4) on a fine, regular mesh. Notice that there is a small overshoot in px at the mould gate, since the boundary condition imposes px ðx ¼ 0:5Þ ¼ 1. Fig. 6 shows the normal pressure gradient at the wall x ¼ 0 where the overshoot at y ¼ 0 can be clearly seen. We also observe an oscillation around the points of discontinuity in qp=qx which is present even in the high resolution solution (N ¼ 337). As a second example we consider the case of Dirichlet boundary conditions at the inlet gate (p ¼ 1). The left part of Fig. 7 shows the RBF solution in the RBF centers, while the right part shows the RBF solution in a fine, regular mesh. Notice the oscillation around p ¼ 1 in the solution at the inlet gate.
qjrpðnþ1Þ jm qpðnþ1Þ qx qx
qjrpðnþ1Þ jm qpðnþ1Þ . qy qy
ð20Þ
The residual of the Neumann boundary conditions is similarly computed. We stop the iteration when maxðjðRðnþ1Þ ðx¯ i ÞÞjÞ; i ¼ 1; . . . ; N p p103 . As a first numerical experiment we consider the case of a Neumann boundary condition at the inlet gate with the following parabolic profile; qp=qx ¼ 102 ðy 0:5Þ2 1. The left of Fig. 3 shows, for different values of the shape parameter, ck ¼ const., how the maximum residual at the collocation centers decreases during the iterative process. It can be observed that the rate of convergence increases with decreasing shape parameter. However, this fact only implies that the nonlinear set of algebraic equations
y
1
0.5
Interior nodes: 50 Boundary nodes: 28 RBF centers: 106
0 0
1
x
867
Fig. 2. Collocation nodes (circles) and RBF centers (crosses).
0.1 ∫ |R(x)| dx
max(abs(R(xk))
100
10-2 c=0.05
c=0.005
10-4
c=0.01 0
10 n
0
20
Fig. 3. Left: convergence of iterative solution procedure for Eq. (19). Right:
0.01
0.05
0.1
0.2
c
R
jRðxÞÞj ¯ d x¯ as a function of c for N ¼ 106 (solid) and N ¼ 337 (dashed).
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
868
0.2
p
R
5
0 0
0
0 y 0.5 0
0.2 x
1
x
0.4 0
1
y
Fig. 4. Left: RBF solution in a fine regular mesh (N ¼ 106; c ¼ 0:04). Right: residuals. Asterisks show the location of RBF centers.
1 1 0 0 0 1 y
0.5 0
x
0 y y
0.5 0
xx
1
Fig. 5. Left: qp=qx in a fine regular mesh (N ¼ 106; c ¼ 0:04). Right: qp=qy.
(Dirichlet boundary condition), the root mean square error is 3:4 103 . 1
4. Front displacement
0
0
1 y
Fig. 6. Imposed pressure gradient at x ¼ 0 (solid), pressure gradient from RBF solution for N ¼ 106 (dashed) and N ¼ 337 (dotted).
In this case the residuals near the inlet gate are even higher, R since there is a discontinuity in first derivatives. In fact, O jRðxÞj dx ¼ 0:54. As a final check of the accuracy of the solution we have computed the RBF solution for the Newtonian case (m ¼ 0), where we can compute the exact solution. In the first example (Neumann boundary condition at the inlet gate), the root mean square error of the solution at the collocation nodes is 5:6 104 . In the second example
Once the pressure field is known, Eq. (6) provides the average velocity. Of particular relevance are the velocities in the advancing front which are needed to update its location. Marker particle techniques [33] and volume-of-fluid techniques [34] are two methods which have been often used to model moving fronts. However, they have severe drawbacks which hamper their use. Since the seminal work of Osher and Sethian [24], level set methods and fast marching methods have been shown to be accurate and numerically efficient techniques, which do not have any of the drawbacks of traditional methods, and which have been successfully applied in a great variety of applications [35,33]. The level set method considers the moving front as the zero level set of a higher dimensional function cðx; tÞ. The evolution of the front is obtained by evolving the level set function with the simple advection equation, ct þ V rc ¼ 0
(21)
which defines the motion of the front cðx; tÞ ¼ 0. It is an Eulerian formulation since the front is captured by the implicit function c. To solve (21) requires the velocity field V, and an initial level set function cðx; 0Þ with the property that its zero level set corresponds to the initial front.
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
869
1
1
p
p
0
0 0
0
1
x 0.5 0
y
x
1 0.5 0
y
Fig. 7. Left: RBF solution at collocation nodes (N ¼ 106). Right: RBF solution in a fine, regular mesh.
Otherwise, the initial level set function is quite arbitrary although a signed distance function is often used. In order to correctly evolve the front, the velocity field V has to be defined throughout the domain or, at least, within a narrow band surrounding it. The front, defined by the zero level set of the solution of (21), will move correctly provided that the velocity of the front is correct. In the present case, the velocity field is only defined in the filled region of the mould and, therefore, it has to be reasonably extended to the whole mould domain. To this end, we compute the velocity of collocation centers located in a narrow band surrounding the front using Eq. (6). Then, extrapolate the velocity field to the rest of the domain. This extrapolation is done using an RBF technique so that the resulting velocities are smooth but not physically correct [31]. However, for the purpose of evolving the front, it is only necessary that the velocity is correct near the front and, therefore, this velocity field is perfectly adequate. Eq. (21) is solved in a uniform mesh covering the whole domain by means of the upwind method, 8 n < ðci;j cni1;j Þ if u¯ i;j X0 nþ1 n ci;j ¼ ci;j nx u¯ ij : ðcniþ1;j cni;j Þ if u¯ i;j o0 8 n < ðci;j cni;j1 Þ if v¯ i;j X0; ny v¯ ij : ðcni;jþ1 cni;j Þ if v¯ i;j o0; where nx ¼ Dt=Dx, ny ¼ Dt=Dy, and where subscripts i (j) refer to the horizontal (vertical) discretization and superscript n to the time discretization (t ¼ n Dt). As boundary conditions we use extrapolation of the level set function. Thus, for instance, at the left boundary (i ¼ 1), cn0;j ¼ 2cn1;j cn2;j and, therefore, n cnþ1 ¯ 1j ðcn2;j cn1;j Þ 1;j ¼ c1;j nx u 8 n < ðc1;j cn1;j1 Þ if v¯ 1;j X0; ny v¯ 1j : ðcn1;jþ1 cn1;j Þ if v¯ 1;j o0:
The time step is defined by a CFL condition Dt ¼
0:9 . maxfj¯uij j=Dx þ j¯vij j=Dyg
(22)
At each time step, the level set function is updated from (21) and a set of points located in the new location of the front are computed by interpolating the level set function in the finite difference mesh. At the next time step, these points become RBF collocation centers, where the boundary condition p ¼ 0 is enforced, while the old boundary centers become interior RBF centers. However, a filtering operation is needed to eliminate points in the front that are too much close to each other or too close to some interior point, and that, therefore, will make the computation of the pressure very ill conditioned. This is carried out by computing the distance of each point to its nearest neighbor, and deleting the point if that distance is below a predefined threshold. In any case, as the front advances, the number of RBF centers increases and therefore the size and condition number of the linear system resulting from Eq. (19) increases. The upper left of Fig. 8 shows the finite difference mesh (30 30, Dx ¼ Dy ¼ 0:6666) and the initial set of RBF centers used. The upper right shows the level set function at t ¼ 0 (signed distance function) together with its evolution at time step t ¼ 0:17. The zero level set of each surface determines the location of the advancing front. The lower left of the figure shows a profile of the initial level set function (solid curve) at x ¼ 0:4 and the same profile at a later time (t ¼ 0:17). The intersection of these curves with the line c ¼ 0 defines the location of the front at each time. The lower right of the figure shows the zero level set (advancing front) at t ¼ 0 and 0:17. The left part of Fig. 9 shows the evolution of the front in a rectangular mould. A good way to assess the accuracy of the method is to check if mass conservation is satisfied. At each time step, the mass (area) injected can be computed by numerically integrating u¯ along the inlet gate using the formulas in the Appendix. Also, the total area filled by fluid can be computed since we know the location of the front at each time step. The right part of Fig. 9 compares the total area injected as a function of time with the area filled by fluid. There is good agreement which shows the accuracy of the method. The results are worse in the case of Dirichlet boundary conditions (p ¼ 1) at the inlet gate. This should be expected since in the case of Neumann boundary condition there are discontinuities in second derivatives which produce large residuals near the inlet
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
870
2
Ψ
y
1
0 0
1 0
1
y
2
0
2
1
0
x
x 2 x=0.4
y
ψ
1
0
0
0
1
0
1
y
2
x
Fig. 8. Upper left: uniform finite difference mesh and initial distribution of RBF centers (crosses) and collocation nodes (circles). Upper right: level set function at t ¼ 0 and 2:02. Lower left: level set function profile cð0:4; y; 0Þ (solid line) and cð0:4; y; 2:02Þ (dashed line). Lower right zero level sets at t ¼ 0 and 2:02; cðx; 0Þ ¼ 0 and cðx; 2:02Þ ¼ 0.
2
y
A
1
1 0
0 0
1 x
2
0
10 t
20
Fig. 9. Left: front location at times 0, 2.02, 4.27, 6.21, 8.30, 10.08, 12.04, 14.24 and 16.24. Right: total area injected (solid line) and total area filled (dashed line) as a function of time.
gate, while in the case of Dirichlet boundary conditions the discontinuities are in first derivatives and therefore much larger residuals occur. To test the suitability of the method, we have used the rectangular mould shown in Fig. 1, but with an elliptic inset placed in the center of the mould. As before, we use a Neumann boundary condition at the inlet gate with the parabolic profile; qp=qx ¼ 102 ðy 0:5Þ2 1. The upper left of Fig. 10 shows the domain, and the initial location of RBF centers and collocation nodes. The other figures show
the evolution of the filling fluid and the RBF centers for t ¼ 6:1; 16:8 and 24:1, respectively. These RBF centers and collocation nodes are used to solve Eq. (13) with the RBF method. As time progresses, the number of RBF centers increases (106 at t ¼ 0, 367 at t ¼ 6:1, 733 at t ¼ 16:8 and 961 at t ¼ 24:1). Thus, the resulting linear system grows and the computing time per step increases. Notice in the lower left of Fig. 10 that the front appears not well defined. This is due to the filtering operation needed to eliminate nodes which are too close to each
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
other, which often results in eliminating nodes that are located in the advancing front. Notice also in the lower right of the figure, how the method is able to handle properly the collision of the two fronts behind the elliptic
871
obstacle. These changes in topology are difficult to handle with other methods. Fig. 11 shows the evolution of the level set function cðx; tÞ which is obtained by solving (21). For t ¼ 0 (upper
1
y
y
1
0
0
0
1 x
2
1
0
1 x
2
0
1 x
2
y
y
1
0
0
0
1 x
2
Fig. 10. Distribution of RBF centers and collocation nodes for t ¼ 0 (upper left), t ¼ 6:1 (upper right), t ¼ 16:8 (lower left) and t ¼ 24:1 (lower right).
1
1
0.8
1
0.8
1
0.6
0.6
Ψ
0.2
0
0.4
Ψ
0.4
0
0.2
0
2
1 1
0
y
0
0
2
1
x
y
1
0 0
x 1
1
0.8
0.8
1
1
0.6
0.6 0.4
0.4
0.2
Ψ
Ψ
0.2
0
0
2
1 y
1
0 0
x
0
0
2
1 y
1
0 0
x
Fig. 11. Level set function evolution: t ¼ 0 (upper left), t ¼ 6:1 (upper right), t ¼ 16:8 (lower left) and t ¼ 24:1 (lower right).
ARTICLE IN PRESS 872
F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
0.3
p
p
0.3
0
0
0
0 1 x
2
0
y
1
x
0.3
1 2
0
2
0
y
1
p
p
0.3
0 0
0 0
1 x
1 2
0
x
y
1
y
1
Fig. 12. RBF solution for pressure at times t ¼ 0 (upper left), t ¼ 6:1 (upper right), t ¼ 16:8 (lower left) and t ¼ 21:9 (lower right).
1
y
left) this function is initialized to the signed distance from the front. Then, this function is advected with the velocity field obtained by RBF interpolation of velocities in the centers located in the narrow band. No re-initialization has been necessary. Notice in the lower right figure that when two fronts collide and produce a change in topology, no difficulty is encountered. The front is always the zero level set of function c, and the two fronts coalesce into a new smooth front. However, the level set function has become very flat and, therefore, small errors in its computation may lead to high errors in the location of the front. This can be avoided by re-initializing the level set function to a signed distance function, when minðrcÞ becomes smaller than a predefined value. Fig. 12 shows the evolution with time of the RBF solution for pressure. The solution is computed from Eq. (15) in all the nodes of a fine, equispaced grid. However, it is only valid in the filled region of the mould which corresponds to nodes where pX0. Notice that the pressure in the inlet gate increases with time, since the pressure gradient is constant, and the distance to the advancing front (where p ¼ 0) increases. When the front reaches the elliptic obstacle, mass conservation implies that the velocities around it have to be larger, and therefore a steepening of the pressure profile is clearly seen. Fig. 13 shows the evolution of the advancing front as time progresses. For clarity, consequent time steps are plotted in different colors (red, blue). It shows the nodes located in the front for different time steps spaced approximately every 3 time units. Although in this case, the velocity of the advancing front is approximately constant, it should be pointed out that as time progresses, the computational cost for advancing a given length increases because the number of RBF centers increases.
0
0
1 x
2
Fig. 13. Front evolution. From left to right: t ¼ 0, 3:258, 6:261, 9:134, 12:08, 15:16, 18:47, 21:10, 24:07, respectively.
5. Conclusions Numerical modeling of injection moulding requires the solution of a moving boundary problem described by a 2D, non-linear elliptic equation for the pressure. We have used the asymmetric RBF collocation method of Kansa, which reduces the problem to the solution of a nonlinear set of algebraic equations. These equations are solved using an iterative procedure based on a simple linearization of the resulting equations. Although the method works well for the problems analyzed in this paper, we are considering an alternative approach based on iteratively solving a linear elliptic equation for the pressure correction.
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
We have found that large errors occur, near the mould entrance, where a discontinuity in the boundary condition exists. However, these large errors are eliminated by using the method proposed by Fedoseyev [32]: in boundary nodes, impose collocation of both the PDE and the boundary condition. We have also found that, near the mould entrance, the residuals in a fine grid can be quite large (specially in the case of Dirichlet boundary conditions which introduce discontinuities in the first derivatives). These large residuals may lead to erroneous volume sources/sinks in the interior which result in violation of mass conservation. In fact, since the method uses the strong form (collocation) approach, conservation is not guaranteed. A possible alternative would be to use the weak form (Galerkin) approach [36,37] to improve mass conservation in all cases. Acknowledgments This work has been supported by the Spanish MEC Grants FLUMRES DPI2002-04550-c07-03 and MAT2005-05730-c02-01. Appendix
fk ðxÞ ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx xk Þ2 þ ðy yk Þ2 þ c2k ,
ð23Þ
qfk x xk ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , qx ðx x Þ2 þ ðy y Þ2 þ c2
ð24Þ
qfk y yk ffi, ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qy 2 ðx x Þ þ ðy y Þ2 þ c2
ð25Þ
k
k
k
k
k
k
q2 f k 1 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qx2 2 ðx xk Þ þ ðy yk Þ2 þ c2k ðx xk Þ2
,
ð26Þ
,
ð27Þ
q2 fk ðx xk Þ ðy yk Þ ¼ , qxqy ððx xk Þ2 þ ðy yk Þ2 þ c2k Þ3=2
(28)
ððx xk Þ2 þ ðy yk Þ2 þ c2k Þ3=2
q2 f k 1 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 qy ðx xk Þ2 þ ðy yk Þ2 þ c2k
ðy yk Þ2 ððx xk Þ2 þ ðy yk Þ2 þ c2k Þ3=2
1 Dfk ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðx xk Þ þ ðy yk Þ2 þ c2k þ
c2k 2
ððx xk Þ þ ðy yk Þ2 þ c2k Þ3=2
,
ð29Þ
Z
qfk dx qx qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ log y yk þ ðx xk Þ2 þ ðy yk Þ2 þ c2k .
873
ð30Þ
References [1] Hardy RL. Multiquadric equations of topography and other irregular surfaces. J Geophys Res 1971;176:1905–15. [2] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics. I. Surface approximations and partial derivative estimates. Comput Math Appl 1990;19:127–45. [3] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics. II. Solutions to parabolic, hyperbolic and elliptic partial differential equations. Comput Math Appl 1990;19:147–61. [4] Madych WR, Nelson SA. Multivariate interpolation and conditionally positive definite functions. Approx Theory Appl 1988;4: 77–89. [5] Madych WR. Miscelaneous error bounds for multiquadric and related interpolators. Comput Math Appl 1992;24(12):121–38. [6] Schaback R. Error estimates and condition numbers for radial basis function interpolation. Adv Comput Math 1995;3:251–64. [7] Beatson RK, Cherrie JB, Mouat CT. Fast fitting of radial basis functions: methods based on preconditioned GMRES iteration. Adv Comput Math 1999;11:253–70. [8] Brown D, Ling L, Kansa EJ, Levesley J. On approximate cardinal preconditioning for solving PDEs with radial basis functions. Eng Anal Bound Elem 2005;29:343–53. [9] Ling L. Radial basis functions in scientific computing. PhD thesis, Simon Fraser University, May 2003. [10] Ling L, Kansa EJ. A least-squares preconditioner for radial basis functions collocation methods. Adv Comput Math 1–25 (2004). [11] Kansa EJ, Hon YC. Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations. Comput Math Appl 2000;39:123–37. [12] Zhou X, Hon YC, Li J. Overlapping domain decomposition method by radial basis functions. Appl Numer Math 2003;44:241–55. [13] Shu C, Ding H, Yeo KS. Local radial basis function-based differential quadrature method and its application to solve twodimensional incompressible Navier–Stokes equations. Comput Methods Appl Mech Eng 2003;192:941–54. [14] Ling L, Hon YC. Improved numerical solver for Kansa’s method based on affine space decomposition. Eng Anal Bound Elem 2005;29:1077–85. [15] Ferreira AJM, Martins PALS, Roque CMC. Solving time-dependent engineering problems with multiquadrics. J Sound Vib 2005;280: 595–610. [16] Ferreira AJM, Roque CMC, Jorge RMN, Kansa EJ. Static deformations and vibration analysis of composite and sandwich plates using a layerwise theory and multiquadrics discretization. Eng Anal Bound Elem 2005;29:1104–14. [17] Hon YC, Cheung KF, Mao XZ, Kansa EJ. A multiquadric solution for the shallow water equation. ASCE J Hydraul Eng 1999;125: 524–33. [18] Hon YC, Lu MW, Xue WM, Zhu YM. Multiquadric method for the numerical solution of biphasic mixture model. Appl Math Comput 1997;88:153–75. [19] Leitao VMA. A meshless method for Kirchoff plate bending problems. J Numer Methods Eng 2001;52:1107–30. [20] Wang JG, Liu GR, Lin P. Numerical analysis of blot’s consolidation process by radial point interpolation method. Int J Solids Struct 2002;39:1557–73.
ARTICLE IN PRESS 874
F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
[21] Fasshauer GE. Newton iteration with multiquadrics for the solution of nonlinear PDEs. Comput Math Appl 2002;43:423–38. [22] Hieber CA, Shen SF. Flow analysis of the non-isothermal twodimensional filling process in injection molding. Isr J Technol 1978;16:248–54. [23] Hieber CA, Shen SF. A finite-element/finite difference simulation of the injection-molding filling process. J Non-Newtonian Fluid Mech 1979;7:1–32. [24] Osher S, Sethian JA. Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations. J Comput Phys 1988;79:12–49. [25] Holm EJ, Langtangen HP. A unified finite element model for the injection molding process. Comput Methods Appl Mech Eng 1999;178:413–29. [26] Kwon TH, Park JB. Finite element analysis modeling of powder injection molding filling process including yield stress and slip phenomena. Polym Eng Sci 1995;35:741–53. [27] Verhoyen O, Dupret F. A simplified method for introducing the Cross viscosity law in the numerical simulation of Hele Shaw flow. J Non-Newtonian Fluid Mech 1998;74:25–46. [28] Franke R. Scattered data interpolation: tests of some method. Math Comput 1982;38:181–200. [29] Wertz J, Kansa EJ, Ling L. The role of the multiquadric shape parameters in solving elliptic partial differential equations. Comput Math Appl 2006;51:1335–48.
[30] Fornberg B, Driscoll TA, Wright G, Charles R. Observation of the behavior of radial basis function approximations near boundaries. Comput Math Appl 2002;43:473–90. [31] Bernal B, Kindelan TA. An RBF meshless method for injection molding modelling, Meshfree Methods for Partial Differential Equations III. In: Griebel M, Scweitzer MA, editors. Lecture notes in computational science and engineering. Springer; 2006. pp. 41–56. [32] Fedoseyev AI, Friedman MJ, Kansa EJ. Improved multiquadric method for elliptic partial differential via PDE collocation on the Boundary. Comput Math Appl 2002;43:439–55. [33] Sethian JA. Level set methods and fast marching methods. Evolving interfaces in computational geometry, fluid mechanics, computer vision and materials science. Cambridge: Cambridge University Press; 1999. [34] Noh W, Woodward PC. A simple line interface calculation. In: van de Vooran AI, Zandberger PJ, editors. Proceedings of the fifth international conference on fluid dynamics. Berlin: Springer; 1976. [35] Osher S, Fedkiw R. Level set methods and dynamic implicit surfaces. Berlin: Springer; 2003. [36] Wedland H. Meshless Galerkin approximation using radial basis functions. Math Comput 1999;68:1521–31. [37] Wedland H. Numerical solution of variational problems by radial basis functions. In: Chui CK, Schumaker LL, editors. Approximation theory IX, vol. II: computational aspects. Vanderbilt University Press; 1999. p. 361–8.
Engineering Analysis with Boundary Elements 31 (2007) 863–874 www.elsevier.com/locate/enganabound
RBF meshless modeling of non-Newtonian Hele–Shaw flow Francisco Bernal, Manuel Kindelan Universidad Carlos III de Madrid, Avda. Universidad 30, 28911 Leganes, Spain Received 31 March 2006; accepted 31 January 2007 Available online 20 April 2007
Abstract In this paper we consider the problem of injecting a non-Newtonian fluid into a thin cavity. Using the Hele–Shaw approximation the problem reduces to a moving boundary problem in which the pressure is described by a 2D, nonlinear, elliptic equation. Mesh-free methods are very well suited for the numerical solution of moving boundary problems since no remeshing is needed at each time step to correctly represent the boundary. Among these methods, we have chosen the asymmetric RBF collocation method (Kansa’s method) to compute the pressure distribution. Kansa’s method is truly mesh-free, and is one of the most frequently used methods due to its accuracy and ease of implementation. Once the pressure is known, the velocity at each point in the moving front is computed and, therefore, the location of the front can be updated. To efficiently model the front motion we use the level set method, which is very accurate, and can handle both front collisions and front break-ups without difficulty. r 2007 Elsevier Ltd. All rights reserved. Keywords: Radial basis function; Generalized multiquadries; Level set method
1. Introduction Recently, very intensive efforts have been devoted to develop mesh-free or element-free methods that eliminate the need of element connectivity in the solution of partial differential equations (PDEs). The motivation is to cut down modeling costs in industrial applications by avoiding the labor intensive step of mesh generation. These methods are particularly attractive in problems with moving interfaces since no remeshing is necessary. In this paper we use the radial basis function (RBF) method which was first used by Hardy [1] for interpolation of scattered data and later by Kansa [2,3] for the solution of PDEs. Kansa’s method looks for an approximate solution in the function space spanned by a finite set of RBFs, by using collocation of the PDE in a set of (scattered) centers. In particular, we use multiquadric (MQ) RBFs which have been shown to have exponential convergence [4,5]. However, the price that one has to pay for this increased accuracy is ill-conditioning of the resulting system of Corresponding author.
E-mail address: [email protected] (M. Kindelan). 0955-7997/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2007.01.008
algebraic equations which have to be solved to determine the coordinates of the approximate solution in the radial functions basis. In fact, there is a trade-off between accuracy and ill-conditioning which Schaback explains by a kind of uncertainty relation [6]. During the last years many procedures have been proposed to improve the conditioning problem; approximate cardinal basis function (ACBF) preconditioning [7–10], nonoverlapping [11] and overlapping [12] domain decomposition methods, differential quadrature methods [13], affine space decomposition [14]. The fact is that, although the RBF method has many attractive properties, it has not been yet extensively applied to solve practical problems. Even more, most of the engineering problems [15–20] where the method has been used are linear ones. Recently, Fasshauer [21] showed that Kansa’s method can be easily adapted for nonlinear elliptic PDEs. In this paper we address the problem of injection molding, which is one of the most important industrial processes for the manufacturing of thin, plastic products. The constant demand for higher quality leads to a growing interest in the analysis and prediction of the product’s physical properties through the use of numerical simulations.
ARTICLE IN PRESS 864
F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
The mathematical model of injection molding is a free boundary problem described by conservation equations for mass, momentum and energy. The fluid is assumed to be non-Newtonian and a power-law equation is used to model its rheology. The conservation equations for this problem, can be dramatically simplified assuming that the mould is thin so that a Hele–Shaw approximation [22,23] can be used. In this case the momentum equation is just a 2D, nonlinear, elliptic equation whose solution yields the pressure distribution in the filled region of the mould. From this pressure field, the velocity distribution can be computed and the location of the advancing front can be updated. The front motion is modeled using a level set method [24], which has proved to be very accurate, and which is able to handle both front collisions and front breakups without difficulty. It is based in considering the moving front as the zero level set of a higher dimensional function which is evolved using the advection equation. For simplicity, the advection equation is numerically solved using the upwind method in a regular equispaced grid. 2. Model equations To model the flow of an injected fluid in a mould, we start from the conservation equations for mass, momentum and energy, qr þ r ðrvÞ ¼ 0, qt qu qu qu qu þu þv þw r qt qx qy qz qp qsx qtxy qtxz þ ¼ þ þ , qx qx qy qz qv qv qv qv þu þv þw r qt qx qy qz qtxy qsy qtyz qp þ ¼ þ þ , qy qx qy qz qw qw qw qw þu þv þw r qt qx qy qz qp qtxz qtyz qsz þ ¼ þ þ , qz qx qy qz qT qT qT qT þu þv þw rC p qt qx qy qz 2 q T q2 T q2 T þ 2 þ 2 þ ZF, ¼k qx2 qy qz
ð2Þ ð3Þ
where the inertia terms in the momentum equations have been neglected in comparison to the viscous terms (creeping flow), and where partial derivatives in the horizontal directions (x; y) have been neglected in comparison to derivatives in the vertical direction (small thickness). Integrating the momentum equations along the vertical direction from the mid-plane, and taking into account that by symmetry ðqu=qzÞðz ¼ 0Þ ¼ ðqv=qzÞðz ¼ 0Þ ¼ 0, results in qu z qv z ¼ Lx ; ¼ Ly , qz Z qz Z qp qp Lx ; Ly . qx qy A new integration in the vertical direction yields, Z b Z b Z b qu Lx z~ z~ d~z ¼ d~z ) u ¼ Lx d~z; q~ z Z Z z z z Z b Z b Z b Ly z~ qv z~ d~z ¼ d~z; d~z ) v ¼ Ly q~ z Z Z z z z
ð4Þ ð5Þ
where bðx; yÞ is the half-height and where the no slip boundary condition, uðz ¼ bÞ ¼ vðz ¼ bÞ ¼ 0, has been used. Integrating these equations along the vertical direction results in the following averaged horizontal velocities: Z Z Ly 1 b Lx 1 b u¯ ¼ S; v¯ ¼ S, (6) u dz ¼ v dz ¼ b 0 b 0 b b where Z S¼ 0
b 2
z dz. Z
Thus, the velocity vector is in the direction of rp. Analogously, the continuity equation (1) is integrated along the half-height, Z b qu qv qðb¯uÞ qðb¯vÞ þ þ . (7) dz ¼ qx qy qx qy 0
where F is the dissipation function. Under the assumptions of incompressibility, creeping flow (ðr1 U 1 L=Z1 Þ ðh=LÞ2 51) and small thickness (h=L51), the above equations simplify to, qu qv þ ¼ 0, qx qy
qp q qu qp q qv qp ¼ Z ¼ Z ¼ 0, ; ; qx qz qz qy qz qz qz qT qT qT q2 T rC p þu þv ¼ k 2 þ ZF, qt qx qy qz
ð1Þ
Introducing the values of the average horizontal velocities (6) into the integral form of the continuity equation (7) produces, q qp q qp S S þ ¼ 0. (8) qx qx qy qy This is an elliptic, nonlinear equation whose solution yields the pressure field and from (4)–(5) the velocity field. Notice that, in general, the viscosity is a function of pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi shear stress (_g ¼ u2z þ v2z ) and temperature (Z ¼ ZðT; g_ Þ).
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
From (4)–(5) g_ ¼
Lz Z
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where L jrpj ¼ L2x þ L2y .
(9)
For instance, if Z is given by the power-law model, namely Z ¼ Z0 gðTÞ_gn1
(11)
In the general case, Z depends on temperature and, therefore, the momentum equation (8) is coupled to the energy equation (3). However, in what follows we will neglect the temperature dependence thus decoupling energy and momentum equations. Taking gðTÞ ¼ 1 the flowability is simply given by nbð2 nþ1Þ=n 1=n Z0 ð2 n
.
(12)
þ 1Þ
We will also assume that the mould width is constant (bðx; yÞ ¼ const.), so that Eq. (8) simplifies to, q qp q qp jrpjm jrpjm þ ¼ 0, (13) qx qx qy qy where m ¼ ð1 nÞ=n. The case m ¼ 0 corresponds to a Newtonian fluid and the case ma0 to a non-Newtonian fluid. From (6) the averaged horizontal velocities are given by u¯ ¼ jrpjm
qp ; qx
v¯ ¼ jrpjm
qp . qy
(14)
Eq. (13) has to be solved with the following boundary conditions:
Many different RBF basis functions have been used in the past. Among them, we have chosen the MQ basis functions, fk ðxÞ ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kx xk k2 þ c2k
inlet pressure (or flow rate) at the mould gate, p ¼ 0 in advancing front, normal velocity zero at the walls of the mould, qp=qn ¼ 0.
The solution of Eq. (13) is used in (14) to compute the average velocity field, which is in turn used to advance the front. This model is the Hele–Shaw approximation for injection molding which was introduced by Hieber and Shen [22,23] and later used by many authors (see for instance [25–27]).
first proposed by Hardy [1] and which have been used extensively both for interpolation and for the solution of PDEs. In a thorough study of interpolation methods for 2D scattered data, Franke [28] found that MQ generally performs best. To specify the MQ basis (16), one must specify both the centers xk and the shape parameters ck . The problem of choosing the optimal shape parameters has been investigated by several authors, although there is no yet a widely accepted theory for choosing them. Recently, Wertz et al. [29] reviewed those efforts and carried out many numerical experiments to study the role of generalized MQ shape parameters in the solution of elliptic PDEs. They conclude that the shape parameters for RBF centers located in the boundary should be much greater than the shape parameters in interior RBF centers. Having specified the MQ basis functions, the problem reduces to determining the unknown coefficients ak ðk ¼ 1; . . . ; NÞ. They are computed by collocation in N ¯ ¼ fx¯ 1 ; x¯ 2 ; . . . ; x¯ N g in the domain, O, collocation nodes, X and in its boundary, qO, where either the PDE (13) or the appropriate boundary conditions are imposed. Usually, the collocation nodes, x¯ i and the RBF centers xi coincide, although it is possible to use different collocation nodes and RBF centers (see for instance [30]). For convenience, Eq. (13) can be rewritten as jrpjm Dp þ
qjrpjm qp qjrpjm qp þ ¼ 0, qx qx qy qy
3. Pressure field To solve Eq. (13) we use an RBF based method. The idea is to find an approximate solution in the function space S ¼ spanffðk xk kÞ : xk 2 Xg, where f is a globally supported RBF and X a set on N RBF centers. Thus, N X k¼1
ak fk ðxÞ;
fk ðxÞ fðkx xk kÞ.
(15)
(17)
which for Newtonian fluids (m ¼ 0) reduces to the Laplace equation. Substituting (15) into (17) and using collocation at node xi , results in, jrpi jm Dpi þ
qjrpi jm qpi qjrpi jm qpi þ ¼ 0, qx qx qy qy
where pi ¼
N X
ak fk ðx¯ i Þ,
k¼1
pðxÞ ¼
(16)
(10)
the flowability S becomes [23] Z 1 L 1=n b z 1=n z dz. S¼ L Z0 g 0
S ¼ Lð1nÞ=n
865
N N qpi X qfk qpi X qfk ¼ ðx¯ i Þ; ¼ ðx¯ i Þ, ak ak qx qx qy qy k¼1 k¼1 2 N X q fk q2 f k ak ð x Þ þ ð x Þ , Dpi ¼ ¯ ¯ i i qx2 qy2 k¼1 2 !2 !2 3m=2 N N X X qf qf ik ik 5 , jrpi jm ¼ 4 ak þ ak qx qy k¼1 k¼1
(18)
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
866
2 !2 !2 3m=21 N N X X qjrpi jm qf qf ik ik 5 ¼ m4 ak þ ak qx qx qy k¼1 k¼1 " ! ! N N X X qfik q2 fik ak ak qx qx2 k¼1 k¼1 ! !# N N X X qfik q2 fik þ ak ak , qy qyqx k¼1 k¼1 2 !2 !2 3m=21 N N X X qjrpi jm qf qf ik ik 5 ¼ m4 ak þ ak qy qx qy k¼1 k¼1 " ! ! N N X X qfik q2 fik ak ak qx qxqy k¼1 k¼1 !# ! N N X X qfik q2 fik þ ak ak , qy qy2 k¼1 k¼1 fi;k ¼ fk ðx¯ i Þ;
qfik qfk ðxÞ ðx¯ i Þ; ¼ qx qx
q2 fik q2 fk ðxÞ ¼ ðx¯ i Þ; qx2 qx2
qfik qfk ðxÞ ðx¯ i Þ, ¼ qy qy
q2 fik q2 fk ðxÞ ¼ ðx¯ i Þ, qy2 qy2
q2 fik q2 fk ðxÞ ðx¯ i Þ, ¼ qxqy qxqy (see the Appendix). Thus, Eq. (18) is a nonlinear algebraic equation on the N unknowns ak , which is imposed in N p PDE collocation nodes x¯ i ; i ¼ 1; . . . ; N p . In addition, there are N b algebraic equations resulting from collocation of the boundary conditions on N b boundary nodes x¯ i ; i ¼ 1; . . . ; N b . If N p þ N b ¼ N, the problem reduces to the solution of N nonlinear equations on the N unknowns ak ; k ¼ 1; . . . ; N. 3.1. Numerical method To solve Eq. (18), we use an iterative procedure based on linearizing equation (18) according to, ðnþ1Þ m jrpðnÞ þ i j Dpi
ðnþ1Þ ðnþ1Þ m m qjrpðnÞ qjrpðnÞ i j qpi i j qpi þ ¼ 0, qx qx qy qy (19)
where the superscriptðnÞ refers to the iteration number. Eq. (19) represents a system of N p linear equations on the N unknowns ak ; k ¼ 1; . . . ; N. A similar procedure is used to linearize the equations resulting from Neumann boundary conditions and, therefore, they provide an additional set of N b linear equations. Thus, at each iteration a system of N linear equations has to be solved to compute the N RBF coefficients ak . We consider the rectangular mould shown in Fig. 1, where the injected fluid is shown in dark. At the inlet gate we specify either the pressure (p ¼ 1, Dirichlet) or the normal derivative (qp=qx ¼ f ðyÞ, Neumann). In the moving
p=0
p=1 ∂ p/∂ x=f(y)
∂ p/∂ x=0
Fig. 1. Problem geometry.
front we take p ¼ 0 (Dirichlet) and at the walls we assume that the normal velocity is zero (qp=qn ¼ 0). For the fluid we define m ¼ 0:6 which corresponds to a non-Newtonian fluid with power-law exponent n ¼ 0:625. In a previous paper [31] dealing with the Hele–Shaw problem in the case of Newtonian fluids, we analyzed this problem using PDE collocation at interior nodes and boundary condition collocation at boundary nodes. We found that large errors occur near the mould entrance, where there is a discontinuity in the boundary condition. These inaccuracies near the boundaries is a well-known feature [30] in all RBF approximations. An alternative to improve the accuracy of the method was proposed by Fedoseyev et al. [32]. It is based on the observation that the residual error is typical largest near the boundaries. Thus, they proposed to enforce collocation of the PDE in boundary nodes so that both the boundary condition and the PDE are imposed in those nodes and the residuals dramatically decrease. However, since the number of equations increases, it is necessary to introduce additional RBF centers to match the number of unknown coefficients ak . In [31] we showed that using this method, or the Not-a-Knot method proposed by Fornberg [30], results in great improvements in accuracy. Fig. 2 shows the initial distribution of collocation nodes and RBF centers. There are N o ¼ 50 interior collocation nodes (red circles), N b ¼ 28 boundary collocation nodes (blue and black circles), and N ¼ 106 RBF centers (crosses). The interior and boundary nodes are used both as collocation nodes and as RBF centers. However, for each collocation node in the boundary, a new RBF center is introduced which is located in the direction normal to the boundary at a distance equal to the average distance between centers. PDE collocation is enforced in N p ¼ N o þ N b ¼ 78 nodes, while boundary condition collocation is enforced in N b ¼ 28 nodes. Thus, there is a total of N ¼ 106 nonlinear equations for the N ¼ 106 unknowns ak ; k ¼ 1; . . . ; N. The N p PDE collocation nodes are taken from a finite element mesh.
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
To start the iterations we use as a first approximation the solution of the Laplace equation (m ¼ 0). Then we proceed with the iterations checking the convergence of the procedure by computing the residual at the collocation nodes, Rðnþ1Þ ðx¯ i Þ ¼ jrpðnþ1Þ jm Dpðnþ1Þ þ þ
resulting from RBF collocation is solved faster than the lower shape parameter of the RBF basis functions. But does this mean that the error of the solution decreases with decreasing shape parameter? To check this, we use Eq. (20) to compute the residual, not at the collocation nodes, x¯ i , but at the nodes of a fine, equidistant mesh. If the computed RBF solution is correct, the residual should be zero at Rany point. To compare the accuracy of the solution we use O jRðxÞj dx which is shown in the right of Fig. 3 as a function of the shape parameter. Notice that the residual decreases with increasing shape parameter until a minimum is reached at approximately c ¼ 0:04. If the number of RBF centers increases, the value of the optimal shape parameter decreases and the residual decreases. This is shown by the dashed line which corresponds to N ¼ 337 RBF centers. A similar behavior was observed in the linear case (m ¼ 0, see [31]). The left part of Fig. 4 shows the RBF solution (15) on a fine, regular mesh. The right part shows the corresponding residuals and the location of the RBF centers. Notice that the residuals are different from zero only in the entrance region, where the discontinuity in the second derivatives cannot be correctly reproduced in the functional space spanned by the smooth RBF basis functions. Notice also, that the residuals are zero at the RBF centers, since they are also used as collocation nodes, but they exhibit large oscillations between them. Fig. 5 shows the first derivatives of the RBF solution (4) on a fine, regular mesh. Notice that there is a small overshoot in px at the mould gate, since the boundary condition imposes px ðx ¼ 0:5Þ ¼ 1. Fig. 6 shows the normal pressure gradient at the wall x ¼ 0 where the overshoot at y ¼ 0 can be clearly seen. We also observe an oscillation around the points of discontinuity in qp=qx which is present even in the high resolution solution (N ¼ 337). As a second example we consider the case of Dirichlet boundary conditions at the inlet gate (p ¼ 1). The left part of Fig. 7 shows the RBF solution in the RBF centers, while the right part shows the RBF solution in a fine, regular mesh. Notice the oscillation around p ¼ 1 in the solution at the inlet gate.
qjrpðnþ1Þ jm qpðnþ1Þ qx qx
qjrpðnþ1Þ jm qpðnþ1Þ . qy qy
ð20Þ
The residual of the Neumann boundary conditions is similarly computed. We stop the iteration when maxðjðRðnþ1Þ ðx¯ i ÞÞjÞ; i ¼ 1; . . . ; N p p103 . As a first numerical experiment we consider the case of a Neumann boundary condition at the inlet gate with the following parabolic profile; qp=qx ¼ 102 ðy 0:5Þ2 1. The left of Fig. 3 shows, for different values of the shape parameter, ck ¼ const., how the maximum residual at the collocation centers decreases during the iterative process. It can be observed that the rate of convergence increases with decreasing shape parameter. However, this fact only implies that the nonlinear set of algebraic equations
y
1
0.5
Interior nodes: 50 Boundary nodes: 28 RBF centers: 106
0 0
1
x
867
Fig. 2. Collocation nodes (circles) and RBF centers (crosses).
0.1 ∫ |R(x)| dx
max(abs(R(xk))
100
10-2 c=0.05
c=0.005
10-4
c=0.01 0
10 n
0
20
Fig. 3. Left: convergence of iterative solution procedure for Eq. (19). Right:
0.01
0.05
0.1
0.2
c
R
jRðxÞÞj ¯ d x¯ as a function of c for N ¼ 106 (solid) and N ¼ 337 (dashed).
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
868
0.2
p
R
5
0 0
0
0 y 0.5 0
0.2 x
1
x
0.4 0
1
y
Fig. 4. Left: RBF solution in a fine regular mesh (N ¼ 106; c ¼ 0:04). Right: residuals. Asterisks show the location of RBF centers.
1 1 0 0 0 1 y
0.5 0
x
0 y y
0.5 0
xx
1
Fig. 5. Left: qp=qx in a fine regular mesh (N ¼ 106; c ¼ 0:04). Right: qp=qy.
(Dirichlet boundary condition), the root mean square error is 3:4 103 . 1
4. Front displacement
0
0
1 y
Fig. 6. Imposed pressure gradient at x ¼ 0 (solid), pressure gradient from RBF solution for N ¼ 106 (dashed) and N ¼ 337 (dotted).
In this case the residuals near the inlet gate are even higher, R since there is a discontinuity in first derivatives. In fact, O jRðxÞj dx ¼ 0:54. As a final check of the accuracy of the solution we have computed the RBF solution for the Newtonian case (m ¼ 0), where we can compute the exact solution. In the first example (Neumann boundary condition at the inlet gate), the root mean square error of the solution at the collocation nodes is 5:6 104 . In the second example
Once the pressure field is known, Eq. (6) provides the average velocity. Of particular relevance are the velocities in the advancing front which are needed to update its location. Marker particle techniques [33] and volume-of-fluid techniques [34] are two methods which have been often used to model moving fronts. However, they have severe drawbacks which hamper their use. Since the seminal work of Osher and Sethian [24], level set methods and fast marching methods have been shown to be accurate and numerically efficient techniques, which do not have any of the drawbacks of traditional methods, and which have been successfully applied in a great variety of applications [35,33]. The level set method considers the moving front as the zero level set of a higher dimensional function cðx; tÞ. The evolution of the front is obtained by evolving the level set function with the simple advection equation, ct þ V rc ¼ 0
(21)
which defines the motion of the front cðx; tÞ ¼ 0. It is an Eulerian formulation since the front is captured by the implicit function c. To solve (21) requires the velocity field V, and an initial level set function cðx; 0Þ with the property that its zero level set corresponds to the initial front.
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
869
1
1
p
p
0
0 0
0
1
x 0.5 0
y
x
1 0.5 0
y
Fig. 7. Left: RBF solution at collocation nodes (N ¼ 106). Right: RBF solution in a fine, regular mesh.
Otherwise, the initial level set function is quite arbitrary although a signed distance function is often used. In order to correctly evolve the front, the velocity field V has to be defined throughout the domain or, at least, within a narrow band surrounding it. The front, defined by the zero level set of the solution of (21), will move correctly provided that the velocity of the front is correct. In the present case, the velocity field is only defined in the filled region of the mould and, therefore, it has to be reasonably extended to the whole mould domain. To this end, we compute the velocity of collocation centers located in a narrow band surrounding the front using Eq. (6). Then, extrapolate the velocity field to the rest of the domain. This extrapolation is done using an RBF technique so that the resulting velocities are smooth but not physically correct [31]. However, for the purpose of evolving the front, it is only necessary that the velocity is correct near the front and, therefore, this velocity field is perfectly adequate. Eq. (21) is solved in a uniform mesh covering the whole domain by means of the upwind method, 8 n < ðci;j cni1;j Þ if u¯ i;j X0 nþ1 n ci;j ¼ ci;j nx u¯ ij : ðcniþ1;j cni;j Þ if u¯ i;j o0 8 n < ðci;j cni;j1 Þ if v¯ i;j X0; ny v¯ ij : ðcni;jþ1 cni;j Þ if v¯ i;j o0; where nx ¼ Dt=Dx, ny ¼ Dt=Dy, and where subscripts i (j) refer to the horizontal (vertical) discretization and superscript n to the time discretization (t ¼ n Dt). As boundary conditions we use extrapolation of the level set function. Thus, for instance, at the left boundary (i ¼ 1), cn0;j ¼ 2cn1;j cn2;j and, therefore, n cnþ1 ¯ 1j ðcn2;j cn1;j Þ 1;j ¼ c1;j nx u 8 n < ðc1;j cn1;j1 Þ if v¯ 1;j X0; ny v¯ 1j : ðcn1;jþ1 cn1;j Þ if v¯ 1;j o0:
The time step is defined by a CFL condition Dt ¼
0:9 . maxfj¯uij j=Dx þ j¯vij j=Dyg
(22)
At each time step, the level set function is updated from (21) and a set of points located in the new location of the front are computed by interpolating the level set function in the finite difference mesh. At the next time step, these points become RBF collocation centers, where the boundary condition p ¼ 0 is enforced, while the old boundary centers become interior RBF centers. However, a filtering operation is needed to eliminate points in the front that are too much close to each other or too close to some interior point, and that, therefore, will make the computation of the pressure very ill conditioned. This is carried out by computing the distance of each point to its nearest neighbor, and deleting the point if that distance is below a predefined threshold. In any case, as the front advances, the number of RBF centers increases and therefore the size and condition number of the linear system resulting from Eq. (19) increases. The upper left of Fig. 8 shows the finite difference mesh (30 30, Dx ¼ Dy ¼ 0:6666) and the initial set of RBF centers used. The upper right shows the level set function at t ¼ 0 (signed distance function) together with its evolution at time step t ¼ 0:17. The zero level set of each surface determines the location of the advancing front. The lower left of the figure shows a profile of the initial level set function (solid curve) at x ¼ 0:4 and the same profile at a later time (t ¼ 0:17). The intersection of these curves with the line c ¼ 0 defines the location of the front at each time. The lower right of the figure shows the zero level set (advancing front) at t ¼ 0 and 0:17. The left part of Fig. 9 shows the evolution of the front in a rectangular mould. A good way to assess the accuracy of the method is to check if mass conservation is satisfied. At each time step, the mass (area) injected can be computed by numerically integrating u¯ along the inlet gate using the formulas in the Appendix. Also, the total area filled by fluid can be computed since we know the location of the front at each time step. The right part of Fig. 9 compares the total area injected as a function of time with the area filled by fluid. There is good agreement which shows the accuracy of the method. The results are worse in the case of Dirichlet boundary conditions (p ¼ 1) at the inlet gate. This should be expected since in the case of Neumann boundary condition there are discontinuities in second derivatives which produce large residuals near the inlet
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
870
2
Ψ
y
1
0 0
1 0
1
y
2
0
2
1
0
x
x 2 x=0.4
y
ψ
1
0
0
0
1
0
1
y
2
x
Fig. 8. Upper left: uniform finite difference mesh and initial distribution of RBF centers (crosses) and collocation nodes (circles). Upper right: level set function at t ¼ 0 and 2:02. Lower left: level set function profile cð0:4; y; 0Þ (solid line) and cð0:4; y; 2:02Þ (dashed line). Lower right zero level sets at t ¼ 0 and 2:02; cðx; 0Þ ¼ 0 and cðx; 2:02Þ ¼ 0.
2
y
A
1
1 0
0 0
1 x
2
0
10 t
20
Fig. 9. Left: front location at times 0, 2.02, 4.27, 6.21, 8.30, 10.08, 12.04, 14.24 and 16.24. Right: total area injected (solid line) and total area filled (dashed line) as a function of time.
gate, while in the case of Dirichlet boundary conditions the discontinuities are in first derivatives and therefore much larger residuals occur. To test the suitability of the method, we have used the rectangular mould shown in Fig. 1, but with an elliptic inset placed in the center of the mould. As before, we use a Neumann boundary condition at the inlet gate with the parabolic profile; qp=qx ¼ 102 ðy 0:5Þ2 1. The upper left of Fig. 10 shows the domain, and the initial location of RBF centers and collocation nodes. The other figures show
the evolution of the filling fluid and the RBF centers for t ¼ 6:1; 16:8 and 24:1, respectively. These RBF centers and collocation nodes are used to solve Eq. (13) with the RBF method. As time progresses, the number of RBF centers increases (106 at t ¼ 0, 367 at t ¼ 6:1, 733 at t ¼ 16:8 and 961 at t ¼ 24:1). Thus, the resulting linear system grows and the computing time per step increases. Notice in the lower left of Fig. 10 that the front appears not well defined. This is due to the filtering operation needed to eliminate nodes which are too close to each
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
other, which often results in eliminating nodes that are located in the advancing front. Notice also in the lower right of the figure, how the method is able to handle properly the collision of the two fronts behind the elliptic
871
obstacle. These changes in topology are difficult to handle with other methods. Fig. 11 shows the evolution of the level set function cðx; tÞ which is obtained by solving (21). For t ¼ 0 (upper
1
y
y
1
0
0
0
1 x
2
1
0
1 x
2
0
1 x
2
y
y
1
0
0
0
1 x
2
Fig. 10. Distribution of RBF centers and collocation nodes for t ¼ 0 (upper left), t ¼ 6:1 (upper right), t ¼ 16:8 (lower left) and t ¼ 24:1 (lower right).
1
1
0.8
1
0.8
1
0.6
0.6
Ψ
0.2
0
0.4
Ψ
0.4
0
0.2
0
2
1 1
0
y
0
0
2
1
x
y
1
0 0
x 1
1
0.8
0.8
1
1
0.6
0.6 0.4
0.4
0.2
Ψ
Ψ
0.2
0
0
2
1 y
1
0 0
x
0
0
2
1 y
1
0 0
x
Fig. 11. Level set function evolution: t ¼ 0 (upper left), t ¼ 6:1 (upper right), t ¼ 16:8 (lower left) and t ¼ 24:1 (lower right).
ARTICLE IN PRESS 872
F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
0.3
p
p
0.3
0
0
0
0 1 x
2
0
y
1
x
0.3
1 2
0
2
0
y
1
p
p
0.3
0 0
0 0
1 x
1 2
0
x
y
1
y
1
Fig. 12. RBF solution for pressure at times t ¼ 0 (upper left), t ¼ 6:1 (upper right), t ¼ 16:8 (lower left) and t ¼ 21:9 (lower right).
1
y
left) this function is initialized to the signed distance from the front. Then, this function is advected with the velocity field obtained by RBF interpolation of velocities in the centers located in the narrow band. No re-initialization has been necessary. Notice in the lower right figure that when two fronts collide and produce a change in topology, no difficulty is encountered. The front is always the zero level set of function c, and the two fronts coalesce into a new smooth front. However, the level set function has become very flat and, therefore, small errors in its computation may lead to high errors in the location of the front. This can be avoided by re-initializing the level set function to a signed distance function, when minðrcÞ becomes smaller than a predefined value. Fig. 12 shows the evolution with time of the RBF solution for pressure. The solution is computed from Eq. (15) in all the nodes of a fine, equispaced grid. However, it is only valid in the filled region of the mould which corresponds to nodes where pX0. Notice that the pressure in the inlet gate increases with time, since the pressure gradient is constant, and the distance to the advancing front (where p ¼ 0) increases. When the front reaches the elliptic obstacle, mass conservation implies that the velocities around it have to be larger, and therefore a steepening of the pressure profile is clearly seen. Fig. 13 shows the evolution of the advancing front as time progresses. For clarity, consequent time steps are plotted in different colors (red, blue). It shows the nodes located in the front for different time steps spaced approximately every 3 time units. Although in this case, the velocity of the advancing front is approximately constant, it should be pointed out that as time progresses, the computational cost for advancing a given length increases because the number of RBF centers increases.
0
0
1 x
2
Fig. 13. Front evolution. From left to right: t ¼ 0, 3:258, 6:261, 9:134, 12:08, 15:16, 18:47, 21:10, 24:07, respectively.
5. Conclusions Numerical modeling of injection moulding requires the solution of a moving boundary problem described by a 2D, non-linear elliptic equation for the pressure. We have used the asymmetric RBF collocation method of Kansa, which reduces the problem to the solution of a nonlinear set of algebraic equations. These equations are solved using an iterative procedure based on a simple linearization of the resulting equations. Although the method works well for the problems analyzed in this paper, we are considering an alternative approach based on iteratively solving a linear elliptic equation for the pressure correction.
ARTICLE IN PRESS F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
We have found that large errors occur, near the mould entrance, where a discontinuity in the boundary condition exists. However, these large errors are eliminated by using the method proposed by Fedoseyev [32]: in boundary nodes, impose collocation of both the PDE and the boundary condition. We have also found that, near the mould entrance, the residuals in a fine grid can be quite large (specially in the case of Dirichlet boundary conditions which introduce discontinuities in the first derivatives). These large residuals may lead to erroneous volume sources/sinks in the interior which result in violation of mass conservation. In fact, since the method uses the strong form (collocation) approach, conservation is not guaranteed. A possible alternative would be to use the weak form (Galerkin) approach [36,37] to improve mass conservation in all cases. Acknowledgments This work has been supported by the Spanish MEC Grants FLUMRES DPI2002-04550-c07-03 and MAT2005-05730-c02-01. Appendix
fk ðxÞ ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx xk Þ2 þ ðy yk Þ2 þ c2k ,
ð23Þ
qfk x xk ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , qx ðx x Þ2 þ ðy y Þ2 þ c2
ð24Þ
qfk y yk ffi, ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qy 2 ðx x Þ þ ðy y Þ2 þ c2
ð25Þ
k
k
k
k
k
k
q2 f k 1 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qx2 2 ðx xk Þ þ ðy yk Þ2 þ c2k ðx xk Þ2
,
ð26Þ
,
ð27Þ
q2 fk ðx xk Þ ðy yk Þ ¼ , qxqy ððx xk Þ2 þ ðy yk Þ2 þ c2k Þ3=2
(28)
ððx xk Þ2 þ ðy yk Þ2 þ c2k Þ3=2
q2 f k 1 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 qy ðx xk Þ2 þ ðy yk Þ2 þ c2k
ðy yk Þ2 ððx xk Þ2 þ ðy yk Þ2 þ c2k Þ3=2
1 Dfk ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðx xk Þ þ ðy yk Þ2 þ c2k þ
c2k 2
ððx xk Þ þ ðy yk Þ2 þ c2k Þ3=2
,
ð29Þ
Z
qfk dx qx qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ log y yk þ ðx xk Þ2 þ ðy yk Þ2 þ c2k .
873
ð30Þ
References [1] Hardy RL. Multiquadric equations of topography and other irregular surfaces. J Geophys Res 1971;176:1905–15. [2] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics. I. Surface approximations and partial derivative estimates. Comput Math Appl 1990;19:127–45. [3] Kansa EJ. Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics. II. Solutions to parabolic, hyperbolic and elliptic partial differential equations. Comput Math Appl 1990;19:147–61. [4] Madych WR, Nelson SA. Multivariate interpolation and conditionally positive definite functions. Approx Theory Appl 1988;4: 77–89. [5] Madych WR. Miscelaneous error bounds for multiquadric and related interpolators. Comput Math Appl 1992;24(12):121–38. [6] Schaback R. Error estimates and condition numbers for radial basis function interpolation. Adv Comput Math 1995;3:251–64. [7] Beatson RK, Cherrie JB, Mouat CT. Fast fitting of radial basis functions: methods based on preconditioned GMRES iteration. Adv Comput Math 1999;11:253–70. [8] Brown D, Ling L, Kansa EJ, Levesley J. On approximate cardinal preconditioning for solving PDEs with radial basis functions. Eng Anal Bound Elem 2005;29:343–53. [9] Ling L. Radial basis functions in scientific computing. PhD thesis, Simon Fraser University, May 2003. [10] Ling L, Kansa EJ. A least-squares preconditioner for radial basis functions collocation methods. Adv Comput Math 1–25 (2004). [11] Kansa EJ, Hon YC. Circumventing the ill-conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations. Comput Math Appl 2000;39:123–37. [12] Zhou X, Hon YC, Li J. Overlapping domain decomposition method by radial basis functions. Appl Numer Math 2003;44:241–55. [13] Shu C, Ding H, Yeo KS. Local radial basis function-based differential quadrature method and its application to solve twodimensional incompressible Navier–Stokes equations. Comput Methods Appl Mech Eng 2003;192:941–54. [14] Ling L, Hon YC. Improved numerical solver for Kansa’s method based on affine space decomposition. Eng Anal Bound Elem 2005;29:1077–85. [15] Ferreira AJM, Martins PALS, Roque CMC. Solving time-dependent engineering problems with multiquadrics. J Sound Vib 2005;280: 595–610. [16] Ferreira AJM, Roque CMC, Jorge RMN, Kansa EJ. Static deformations and vibration analysis of composite and sandwich plates using a layerwise theory and multiquadrics discretization. Eng Anal Bound Elem 2005;29:1104–14. [17] Hon YC, Cheung KF, Mao XZ, Kansa EJ. A multiquadric solution for the shallow water equation. ASCE J Hydraul Eng 1999;125: 524–33. [18] Hon YC, Lu MW, Xue WM, Zhu YM. Multiquadric method for the numerical solution of biphasic mixture model. Appl Math Comput 1997;88:153–75. [19] Leitao VMA. A meshless method for Kirchoff plate bending problems. J Numer Methods Eng 2001;52:1107–30. [20] Wang JG, Liu GR, Lin P. Numerical analysis of blot’s consolidation process by radial point interpolation method. Int J Solids Struct 2002;39:1557–73.
ARTICLE IN PRESS 874
F. Bernal, M. Kindelan / Engineering Analysis with Boundary Elements 31 (2007) 863–874
[21] Fasshauer GE. Newton iteration with multiquadrics for the solution of nonlinear PDEs. Comput Math Appl 2002;43:423–38. [22] Hieber CA, Shen SF. Flow analysis of the non-isothermal twodimensional filling process in injection molding. Isr J Technol 1978;16:248–54. [23] Hieber CA, Shen SF. A finite-element/finite difference simulation of the injection-molding filling process. J Non-Newtonian Fluid Mech 1979;7:1–32. [24] Osher S, Sethian JA. Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations. J Comput Phys 1988;79:12–49. [25] Holm EJ, Langtangen HP. A unified finite element model for the injection molding process. Comput Methods Appl Mech Eng 1999;178:413–29. [26] Kwon TH, Park JB. Finite element analysis modeling of powder injection molding filling process including yield stress and slip phenomena. Polym Eng Sci 1995;35:741–53. [27] Verhoyen O, Dupret F. A simplified method for introducing the Cross viscosity law in the numerical simulation of Hele Shaw flow. J Non-Newtonian Fluid Mech 1998;74:25–46. [28] Franke R. Scattered data interpolation: tests of some method. Math Comput 1982;38:181–200. [29] Wertz J, Kansa EJ, Ling L. The role of the multiquadric shape parameters in solving elliptic partial differential equations. Comput Math Appl 2006;51:1335–48.
[30] Fornberg B, Driscoll TA, Wright G, Charles R. Observation of the behavior of radial basis function approximations near boundaries. Comput Math Appl 2002;43:473–90. [31] Bernal B, Kindelan TA. An RBF meshless method for injection molding modelling, Meshfree Methods for Partial Differential Equations III. In: Griebel M, Scweitzer MA, editors. Lecture notes in computational science and engineering. Springer; 2006. pp. 41–56. [32] Fedoseyev AI, Friedman MJ, Kansa EJ. Improved multiquadric method for elliptic partial differential via PDE collocation on the Boundary. Comput Math Appl 2002;43:439–55. [33] Sethian JA. Level set methods and fast marching methods. Evolving interfaces in computational geometry, fluid mechanics, computer vision and materials science. Cambridge: Cambridge University Press; 1999. [34] Noh W, Woodward PC. A simple line interface calculation. In: van de Vooran AI, Zandberger PJ, editors. Proceedings of the fifth international conference on fluid dynamics. Berlin: Springer; 1976. [35] Osher S, Fedkiw R. Level set methods and dynamic implicit surfaces. Berlin: Springer; 2003. [36] Wedland H. Meshless Galerkin approximation using radial basis functions. Math Comput 1999;68:1521–31. [37] Wedland H. Numerical solution of variational problems by radial basis functions. In: Chui CK, Schumaker LL, editors. Approximation theory IX, vol. II: computational aspects. Vanderbilt University Press; 1999. p. 361–8.