Ahsan A. Two Phase Flow, Phase Change and Numerical Modeling
Подождите немного. Документ загружается.
Modeling the Physical Phenomena Involved by Laser Beam – Substance Interaction
9
22
n
2
TL
v
nk
∂ρ
=−
∂
(24)
where:
11
2
k[Wm K ]
−−
⋅⋅ is the heat conductivity that belongs to liquid state,
2
L[J/k
g
] – the
latent melting heat,
3
2
[k
g
/m ]ρ – the mass density that belongs to liquid state, and
n
v[m/s] is the movement speed of the boundary surface,
l
C(t), in the direction of its
external normal vector
n .
The boundary at the t moment is supposed as known, respectively the points
(r(t), 0) and
(0, z(t)) on it. It is enough to determine the points
()
()
rt t,0+Δ and
()
()
0, z t t+Δ in order
to find
l
C(t t)+Δ . In the point (r(t), 0) , (24) yields:
22
r
2
TL
v
rk
∂ρ
=−
∂
2
r
22
T
k
v
Lr
∂
=−
ρ
∂
(25)
where:
()
top
Tr r T
T
rr
+Δ −
∂
=
∂Δ
It obtains:
()
()
r
rt t rt v t+Δ = + ⋅Δ
(26)
In
(0, z(t)) point, (24) yields:
22
z
2
TL
v
zk
∂ρ
=−
∂
2
z
22
T
k
v
Lz
∂
=−
ρ
∂
(27)
where:
()
top
Tz z T
T
zz
+Δ −
∂
=
∂Δ
. It results:
()
()
z
zt t zt v t+Δ = + ⋅Δ
(28)
The new boundary parameters,
(t t)α+Δ
and
(t t)
β
+Δ
, are returned by (26) and (28):
()()
ttrttα+Δ= +Δ,
()()
ttzttβ+Δ= +Δ
(29)
The moment
top
t is the first time when the above procedure is applied.
top
z(t ) 0= and
top
r(t ) 0= at this moment of time. Because the temperature gradient (having the z direction,)
is known in z 0= and r 0= :
()
L
l
z0
T
1
0,0,t
zk
=
∂
=− ϕ
∂
(30)
in (28) results:
()
()
L
top
22
0,0,t
zt t t
L
ϕ
+Δ = Δ
ρ
(31)
Two Phase Flow, Phase Change and Numerical Modeling
10
where
L
ϕ
is the power flow on the processed surface corresponding to the liquid state. In
these conditions, (26) becomes:
()
()
top
top 2
22
TTr
rt t k t
rL
−Δ
+Δ = ⋅Δ
Δ⋅ρ⋅
(32)
The same procedure is applied to find the
v
C(t) boundary, taking into account the latent
heat of vaporization
3
L[J/k
g
] , the mass density corresponding to vapor state
3
3
[k
g
/m ]ρ
and respectively the heat conductivity corresponding to vapor state
11
3
k[Wm K ]
−−
⋅⋅ .
2.4 Digitization of heat equation, boundary and initial conditions
The first step of the mathematical approach is to make the equations dimensionless
(Mazumder, 1991; Pearsica et al., 2008a, 2008c). In heat equation case it will be achieved by
considering the following (
r
∞
and
z
∞
are the studied domain boundaries, where the
material temperature is always equal to the ambient one):
2
a
1
r
r xr, z yr, T Tu, t
K
∞
∞∞
== ==τ
(33)
The heat equation (12) in the new variables
x
,
y
,
τ
, and u yields:
22
1
22
i
uu
1uuK
xx x y K
∂∂
∂∂
++=
∂∂ ∂ ∂τ
,
()
[
]
[
]
x,y 0,1 0,1∈×, 0τ≥ , and
i 1,2,3=
(34)
The initial and limit conditions for the unknown function, u yield:
-
phase 1, for
top
0tt≤<
u(x,
y
,0) 1= , (x,
y
)[0,1][0,1]∈×
(35)
u(1,
y
,) 1τ= ,
y
[0, 1]∈ ,
top
[0, ]τ∈ τ ,
1
top top
2
K
t
r
∞
τ=
(36)
u(x,1, ) 1τ= , x[0,1]∈ ,
top
[0, ]τ∈ τ
(37)
-
phase 2, for
top vap
ttt≤<
top top
u(0,0, ) uτ=
(38)
top 1 top
u(x,
y
,)u(x,
y
,)τ= τ,
(x,
y
)(0,1](0,1]∈×
(39)
where:
top
τ is the τ value when
top
uu,=
top top a
uT/T,= and
1top
u(x,
y
,)τ is the heat
equation solution in according to phase 1.
If
top vap
[, )τ∈ τ τ both solid and liquid phases coexist in material, occupying
s
D()τ and
l
D( )τ
domains respectively, which are separated by a time varying boundary,
l
C( )τ , so
top
u(x,
y
,) uτ= on it. The projection of the domain
l
D( )τ on
y
0= plane is the set
1
{x /x x }≤ . For x 1= and
y
1= respectively, the conditions are:
Modeling the Physical Phenomena Involved by Laser Beam – Substance Interaction
11
u(1,
y
,) u(x,1,) 1τ= τ=
(40)
Phase 2 is going on while
top vap
[, )τ∈ τ τ , where:
vap 1
vap
2
tK
r
∞
⋅
τ=
.
-
phase 3, for
vap
tt≥
vap vap
u(0,0, ) uτ=
(41)
vap 2 vap
u(x,
y
,)u(x,
y
,)τ= τ,
lvap
(x,
y
)D( )\(0,0)∈τ
(42)
vap 1 vap
u(x,
y
,)u(x,
y
,)τ= τ,
svap
(x,
y
)D( )∈τ
(43)
where
2vap
u(x,
y
,)τ is the heating equation solution from phase 2. In this temporal phase all
the three (solid, liquid and vapor) states coexist in material, occupying the domains:
s
D()τ
,
l
D( )τ and
v
D()τ , separated by mobile boundaries
l
C( )τ and
v
C()τ , on which
vap
u(x,
y
,) uτ= . The projection of the domains
l
D( )τ and
v
D()τ on plane
y
0= are the sets:
21
{x /x [x , x ]}∈
and
2
{x /x [0, x ]}∈
. According to phase 3, the conditions on
y
0=
surface
(Neumann type conditions) are:
a.
21
d
xx
r
∞
≤≤:
()
[]
() (
]
()
Vf rc 2
a3
L21
a2
S1
a1
r
x,
y
,,x0,x
Tk
r
x,0, , x x , x
Tk
u
y
rd
x,0, , x x ,
Tk r
d
0, x ,1
r
∞
∞
∞
∞
∞
−ϕτ−ϕ−ϕ∈
⋅
−ϕτ∈
⋅
∂
=
∂
−ϕτ∈
⋅
∈
(44)
b.
21
d
xx
r
∞
≤≤:
()
[]
()
Vf rc 2
a3
L2
a2
r
x,
y
,,x0,x
Tk
u
rd
x,0, , x x ,
yTk r
d
0, x ,1
r
∞
∞
∞
∞
−ϕτ−ϕ−ϕ∈
⋅
∂
=− ϕ τ ∈
∂⋅
∈
(45)
Two Phase Flow, Phase Change and Numerical Modeling
12
c.
2
d
x
r
∞
> :
()
Vf rc
a3
rd
x,y , , x 0,
Tk r
u
y
d
0, x ,1
r
∞
∞
∞
−ϕτ−ϕ−ϕ∈
⋅
∂
=
∂
∈
(46)
Similar Neumann type conditions are settled for temporal phases 1 and 2, accordingly to
their specific parameters.
For x 1
= , and
y
1= respectively, the conditions are given by (40).
2.5 Digitization of equations on separation boundaries
The speed of time variation of separation boundaries,
n
v
, is given by (47), where n is the
external normal vector of the boundary.
e
n
ee
T
k
v
Ln
∂
=−
ρ
⋅∂
,
e2,3=
(47)
For
y
0= and
f
xx= , it results:
eea
rf
ee ee
Tu
kkT
v(x,0)
Lr Lrx
∞
∂∂
⋅
=− =−
ρ
⋅∂
ρ
⋅⋅∂
, e2,3=
(48)
respectively,
1
rf
dr K dx
v(x,0)
dt r d
∞
==
τ
(49)
It results:
ea
ee 1
u
dx k T
dLKx
∂
⋅
=−
τ
ρ
⋅⋅ ∂
, e2,3=
(50)
The
α parameter of separation boundary at τ+Δτ moment is:
ea
ff f
ee 1
u
dx k T
x( ) x( ) x( )
dLKx
∂
⋅
α= τ+Δτ = τ + Δτ= τ − Δτ
τρ⋅⋅∂
, e2,3=
(51)
where:
kf
kf
u u(x ,0) u(x ,0)
xxx
∂−
≈
∂−
(52)
where
k
x ∈ digitization network. For
top
τ=τ and
vap
τ=τ respectively, (52) yields:
11 top
ftop
1
u(x,0) u
u
,x( )0
xx
−
∂
=τ=
∂
(53)
Modeling the Physical Phenomena Involved by Laser Beam – Substance Interaction
13
21 vap
fvap
1
u(x,0) u
u
,x( )0
xx
−
∂
=τ=
∂
(54)
For x 0
= and
f
yy
= , it results:
eea
zf
ee ee
Tu
kkT
v(0,y)
Lz Lr
y
∞
∂∂
⋅
=− =−
ρ
⋅∂
ρ
⋅⋅∂
, e2,3=
(55)
respectively,
1
zf
d
y
dz K
v(0,y)
dt r d
∞
==
τ
(56)
It results:
ea
ee 1
d
y
u
kT
dLK
y
∂
⋅
=−
τ
ρ
⋅⋅ ∂
, e2,3=
(57)
The
β
parameter of separation boundary at
τ+Δτ
moment is:
ea
ff f
ee 1
dy u
kT
y( ) y() y( )
dLKy
∂
⋅
β
=τ+Δτ=τ+Δτ=τ− Δτ
τρ⋅⋅∂
,
e2,3=
(58)
where:
kf
kf
uu(0,
y
)u(0,
y
)
yyy
∂−
≈
∂−
(59)
For
top
τ=τ it results:
Ltopftop
a2
u
r
(0,0, ),
y
()0
yTk
∞
∂
=−
ϕ
ττ=
∂⋅
(60)
For
vap
τ=τ it results:
Vtoprcfvap
a3
u
r
(0,0, ) ,
y
()0
yTk
∞
∂
=−
ϕ
τ−
ϕ
−
ϕ
τ=
∂⋅
(61)
3. Determination of temperature distribution in material
Using the finite differences method, the domain [0, 1] [0, 1]× is digitized by sets of
equidistant points on Ox and Oy directions (Pearsica et al., 2008a, 2008b).
3.1 Digitization of mathematical model equations
In the network points, the partial derivatives will be approximated by:
()
i1,
j
i1,
j
i,j
uu
u
x2x
+−
−
∂
≈
∂Δ
,
()
()
2
i1,
j
i,
j
i1,
j
2
2
i,j
u2uu
u
x
x
+−
−+
∂
≈
∂
Δ
(62)
Two Phase Flow, Phase Change and Numerical Modeling
14
()
i,
j
1i,
j
1
i,j
uu
u
y2y
+−
−
∂
≈
∂Δ
,
()
()
2
i,
j
1i,
j
i,
j
1
2
2
i,j
u2uu
u
y
y
+−
−+
∂
≈
∂
Δ
(63)
()
()
()
i,j i,j
i,j
uu
u
τ+Δτ − τ
∂
≈
∂τ Δτ
(64)
With these approximations, in each inner point of the network the partial derivatives
equations become an algebraic system such:
j,i
1j,i
j
i,i
1j,i
j
i,i
j,1i
j
1i,i
j,i
j
i,i
j,1i
j
1i,i
fucubuauaua =++++
+
−+
+
−
−
(65)
The system coefficients are linear expressions of the partial derivatives equation, computed
in the network points. If there are M and N points on Ox and Oy axis respectively, the
system will include NM
× equations with
()()
1N1M +×+ unknowns. Adding the
conditions for the domain boundaries, the system is determinate.
The implicit method, involving evaluations of the equation terms containing spatial
derivatives at
τΔ+τ moment, is used to obtain the unknown function
()
τ,y,xu distribution
in network points. The option is on this method because there are no restrictions on
choosing the time and spatial steps
()
y,x, ΔΔτΔ . According to this method, an additional
index is introduced, representing the time moment. With these explanations, the heat
equation with finite differences yields:
1n,j,i
e
1
1n,j,i
2
2
2
2
u
K
K
y
u
x
u
x
u
x
1
+
+
τ∂
∂
=
∂
∂
+
∂
∂
+
∂
∂
,
3,2,1e =
(66)
Finally, the algebraic system yields:
()
+
Δ
+λ+
λ+λ+−
Δ
−λ
++++− 1n,j,1i
i
11n,j,i21
e
1
1n,j,1i
i
1
u
x2
x
1u12
K
K
u
x2
x
1
e
1
n,j,i1n,1j,i211n,1j,i21
K
K
uuu ⋅−=λλ+λλ+
+++−
, 3,2,1e =
(67)
where:
()
1
2
x
∂τ
λ=
Δ
,
2
2
x
y
Δ
λ=
Δ
, and
()
io
xx i1 x=+−⋅Δ
. The value
o
x is very close to zero
and it was chosen to avoid the singularity appearing in heating equation at x 0
= . Formally,
this singularity appears because if x 0
= , then
u
0
x
∂
=
∂
too. If x1/MΔ= and
y
1/NΔ=
will result
i1,M1=+ and
j
1,N 1=+.
Equation (67) will be written for
M,2i = and N,2j = . In case 1Nj += and 1Mi += , the
constraints imposed to u are:
1u
1n,1N,i
=
++
, 1u
1n,j,1M
=
++
(68)
Modeling the Physical Phenomena Involved by Laser Beam – Substance Interaction
15
The initial condition is for 0n = :
1u
1,j,i
= ,
()
j,i∀
(69)
For
1i = , 1j ≠ , (67) is still valid, observing that
1n,j,21n,j,0
uu
++
= , because the solution is
symmetrical related to 0x
= . In this case, (67) yields:
()
1
1 2 1,j,n 1 1 2, j,n 1 1 2 1, j 1,n 1 1 2 1,j 1,n 1
e
K
21 u 2u u u
K
++−+++
−+λ+λ +λ +λλ +λλ =
1
1,j,n
e
K
u
K
=− ⋅
, 3,2,1e = (70)
For
1j = , when writing the initial conditions for the boundary 0y = , the temporal phase of
the material must be taken into account. Only the equations corresponding to the third
phase (
vap
tt ≥
) will be presented, because it is the most complex one, all the three states
(solid, liquid and vapor) being taken into account. Similar results were obtained for the
other two phases, in a similar way, accordingly to their influencing parameters. The initial
conditions are:
0vapvapn,1,1
n,uu
0
⋅τΔ=τ=
(71)
()
()
()
0,0\Dy,x,uu
l
ji
2
n,j,i
n,j,i
0
0
∈=
(72)
()
()
sji
1
n,j,i
n,j,i
Dy,x,uu
0
0
∈=
(73)
where
)2(
u is the solution of the problem corresponding to the second temporal phase
(
vaptop
ttt <≤ ).
The boundary corresponding to the
1n
0
+ time moment will be determined hereinafter.
The parameters of the boundary separating the vapor and liquid will be:
()
000
a3
n 1 1,1,n 2,1, n
1v 3
Tk
uu
KLx
+
⋅⋅Δτ
α= −
⋅ρ ⋅ ⋅Δ
,
()
0
n1 V
v3 1
r
0,0,
LK
∞
+
⋅Δτ
β
=
ϕ
τ
ρ⋅ ⋅
(74)
The boundary equation will be:
0
0
0
n1
22
n1
n1
x
y
+
+
+
α
=β−
β
(75)
The boundary separating liquid and solid phases exists at the moment n
0
, as well at a certain
moment n, its parameters being given by (51) and (58). The situation corresponding to the
third phase is illustrated in figure 1.
In the marked points, the heating equation must be changed, because its related partial
derivatives approximation by using finite differences is not possible anymore (the
associated Taylor series in A and B will be used, where
AC a x=Δ and BC b
y
=Δ). In order
Two Phase Flow, Phase Change and Numerical Modeling
16
to know if a point is nearby the boundary previously determined, the points
''
ii
(x ,
y
) and
""
jj
(x , y )
respectively will be obtain by intersecting the network lines
()
io
(x x i 1 x,=+−Δ
()
j
yj
1
y
)=−Δ and the boundary
22
x
y
α
=β−
β
. It results
'
i
y
h1
y
=+
Δ
(floor() + 1) as
solution of the equation
()
'
i
h1
yy
−Δ= , and:
() ()
''
ii
pm
yy
bi,h ,b i,h 1
yy
=− =
ΔΔ
,
() ()()
mpp
bi,h1 1bi,h,bi,h1 1+=− +=
(76)
Fig. 1. The boundaries separating the phases
It results (similar)
"
j
x
k1
x
=+
Δ
as solution of the equation
()
"
j
k1 xx−Δ=, and:
() ()
""
jj
pm
xx
ak,
j
,a k,
j
1
xx
=− =
ΔΔ
,
() ()()
mpp
a k 1,j 1 a k,j , a k 1,j 1+=− +=
(77)
For the regular points (the ones which are not nearby the boundary), the above mentioned
parameters will be:
() () () ()
pmpm
bi,
j
bi,
j
ai,
j
ai,
j
1====
(78)
The unified heating equation in a network point (i, j) will be:
()
pmp
1
11112
i 1,j,n 1 i,j,n 1
mm p i 3 mp i mp mp
1im
12 12 1
fx, n 1 i, j 1 fy, n 1 i , j, n
pm p i mp m pm p 3
ax a a
x
2K22
1u u
a(a a) 2x K aa x aa bb
2x a x
22 K
uuuu
a(a a)x b (b b ) b(b b) K
−+ +
+− +
Δ−
λΔ
λλλλ
−++ ++ −
+
λ+Δ
λλ λλ
−−−=⋅
+++
(79)
Modeling the Physical Phenomena Involved by Laser Beam – Substance Interaction
17
where the coefficients
mpm
a,a,b, and
p
b depends on point (i, j), and
fx i 1,
j
uu
+
= and
f
y
i,
j
1
uu
+
= , if the point where derivatives are approximated is not nearby the boundary. For
i1= , (79) becomes:
pm p
m
1112 1 1
1,j,n 1
epm pm pm i pm mim pi
12 12 1
2, j,n 1 1,j 1,n 1 1,j 1,n 1 1, j,n
pmmp pmp e
aa ax
ax
K2 2 x 2 1
u
Kaa bb aa x aa2axa 2ax
12 2 K
uu uu
ab(bb)b(bb) K
+
+−+++
−Δ
Δ
λλλ λ⋅Δ λ
++ − + −− −
+
λλ λλ
−− − =⋅
++
(80)
There are as well vapor, liquid and solid zones on the boundary
y
0= . Depending on the
position of the intersecting points between boundary and
y
0= , the following situations
may occur:
a.
1
d
x
r
∞
≤ :
()
[]
m
i,1,n 1 fr V i f i 2
a3
ya r
uu x,
y
,,x 0,x
Tk
∞
+
Δ⋅ ⋅
−= ϕ τ ∈
⋅
(81)
()
(
]
i,1,n 1 i,2,n 1 L i i 2 1
a2
yr
uu x,0,,xx,x
Tk
∞
++
Δ⋅
−=ϕτ∈
⋅
(82)
()
i,1,n 1 i,2,n 1 S i i 1
a1
yr
d
uu x,0,,xx,
Tk r
∞
++
∞
Δ⋅
−=ϕτ∈
⋅
(83)
b.
21
d
xx
r
∞
≤<:
()
[]
m
i,1,n 1 fr V i f i 2
a3
ya r
uu x,
y
,,x 0,x
Tk
∞
+
Δ⋅ ⋅
−= ϕ τ ∈
⋅
(84)
()
i,1,n 1 i,2,n 1 L i i 2
a2
yr
d
uu x,0,,xx,
Tk r
∞
++
∞
Δ⋅
−=ϕτ∈
⋅
(85)
c.
2
d
x
r
∞
> :
()
m
i,1,n 1 fr V i f i
a3
ya r
d
uu x,y,,x0,
Tk r
∞
+
∞
Δ⋅ ⋅
−= ϕ τ ∈
⋅
(86)
If
i
d
x
r
∞
> , in all of the three mentioned cases:
i,1,n 1 i,2,n 1
uu
++
= (87)
Two Phase Flow, Phase Change and Numerical Modeling
18
Because the discrete network parameters do not influence the initial moment of laser
interaction with the material, the temperature gradient on z direction was replaced by the
temporal temperature gradient in the initial condition on z 0
= boundary (Draganescu &
Velculescu, 1986). So, the digitized initial condition on z 0
= boundary yields:
()
i,1,n 1 i,1,n i
1
rKd
uu x,0,
kKn
∞
+
⋅τ
=+ ⋅
ϕ
τ
⋅
(88)
The equations system obtained after digitization and boundary determination will be
solved by using an optimized method regarding the solving run time, namely the
column wise method. It is an exact type method, preferable to the direct matrix inversing
method.
3.2 The column wise solving method
From the algebraic system of
(M 1) (N 1)+× +
equations, the minimum dimension will be
chosen as unknowns’ column dimension. It is assumed to be M 1
+ . It is to notice that
writing the system in the point
(i,
j
)
involves as well the points
(i 1,
j
), (i 1,
j
), (i,
j
1),−+ −
and
(i,
j
1)+
(Pearsica et al., 2008a, 2008b). The system and transformed conditions may be
organized, writing in sequence all the equations for each fixed j and variable i, as a vector
system. So, by keeping j constant, it results a relationship between columns j,
j
1−
, and
j
1+
. By denoting [A
j
], [B
j
], and [C
j
] the unknowns coefficients matrixes of the columns j,
j
1−
, and
j
1+
respectively, the system for j constant will be:
jj jj
1
jj
1
j
[A ] {U } [B ] {U } [C ] {U } {F }
−+
⋅+⋅ +⋅ = (89)
where: [X] is a quadratic matrix, {X} is a column vector, {F
j
} is the free terms vector, [A
j
] is a
tri-diagonal matrix whose non-null components are
i,i 1
a
−
,
i,i
a and
i,i 1
a
+
, and [B
j
] and [C
j
]
are diagonal matrixes. The components of matrixes [A
j
], [B
j
] and [C
j
], are the coefficients of
the unified caloric equation written in a point (i,j) of the network, equation (79). The
components of {F
j
} are:
1
i,j i,j,n
e
K
fu
K
=⋅,
j
1≠ , e 1,2,3= (90)
For
j
1= , (90) yields
1
([A ] [I]= - unity matrix, and
11
[B ] [C ] [0]):==
11 1
[A ] {U } {F }⋅= (91)
The components of {F
1
} are computed using the relation:
()
i,1,n i d
1
i,1
i,1,n d
rKd
ux,0,,ii
kKn
f
u,ii
∞
⋅τ
+⋅
ϕ
τ≤
⋅
=
>
(92)
where
d
i is the laser beam limit. Taking into account the relation linking two successive
columns,
j
1
U
−
and
j
U: