Ece - assignment - define the finite volume vi of node i

Question

1. The magnetic field generated from a time-varying electric field in a medium with electric permittivity ε and electric conductivity σ is given by Ampere's Law:

∇ x H¯ = ε ∂E¯/∂t + σE¯.       (1)

If the medium is invariant in z (∂/∂z = 0) and the electric field has only a z component (E_z), then the magnetic field is confined to the 2D x-y plane (i.e., H_z = 0). The above equation is to be discretized using the Finite Volume method on a 2D mesh composed of equilateral triangles of side d as shown below. Suppose the electric field E_z is assigned to each node i in the mesh and time level n, E_z ( x_i , y_i ,t_n ) = Ez_n.

(a) Define the finite volume V_i of node i and indicate the appropriate location assignment of the magnetic field H in the mesh.

(b) By applying Ampere's law to the finite volume V_i of node i, obtain a Finite Volume discrete equation for Eq. (1) that is accurate to 2^nd-order in time. Express all coefficients in the discrete equation in terms of the length d of the triangle and the time step Δt. (Hint: use Stokes' theorem to convert the surface integral of a curl to a contour integral.)

2. A microfluidic channel has a cross-sectional geometry in the x-y plane as shown Figure 2.1. The channel is assumed to have infinite length. Suppose fluid flow in the channel is modeled by the steady-state Navier-Stokes equation for incompressible flow:

(v¯.∇)v¯ - 1/Re∇²v¯ = -∇p,

where v¯ is the flow velocity, p is the pressure and Re is Reynolds number.

(a) If the flow is driven by a constant pressure gradient in the z-direction, ∇p = - p_zzˆ, then the flow velocity can be assumed to have only a z-component, in this case the above Navier-Stokes equation reduces to v¯ = v(x, y)zˆ. Show that in this 

Re ∇ v( x, y) = - pz , (pz = ∂p/∂z) (1)

(b) Show that the Finite Volume discretization of Eq. (1) for a general node i on a triangular mesh can be written as
ni wik ni wik

- v_iΣⁿⁱ_k=1w_ik/l_ik + Σⁿⁱ_k=1w_ik/l_ik = v_k = -Re.p_z A_i , (2)

where v_i = v(x_i, y_i), Ai is the area of the Voronoi polygon associated with node i, ni is the number of neighbour nodes of node i, lik is the length of edge i-k connecting nodes i and k, and wik is the length of the side of the Voronoi polygon perpendicular to edge i-k (see Figure 2.3).

(c) Write a MatLab program to implement the discrete system given by Eq. (2) and solve for the cross-sectional velocity field v(x, y). A triangular mesh for the cross section of the channel is shown in Figure 2.2. The P and T arrays describing the mesh can be downloaded from eClass website (MatLab data file channel_mesh.mat). Assume Re = 1, pz = 1. The boundary condition on the walls of the channel is specified as v = 0 (this is known as the no-slip flow condition). An outline of the MatLab program is given on the next page. Refer to Figure 2.4 for various parameters in the algorithm as well as the formula for calculating the circumcenter of a triangle.

(d) Display the solution v(x, y) on a 3D plot using the trimesh or trisurf function. Hand in a copy of your MatLab program.

 

Outline of the Finite Volume Program

% ---- load mesh files: load channel_mesh
Ne = length(T); Npts = length(P);

% ---- construct matrix A and RHS vector b: A = zeros(Npts,Npts);
b = zeros(Npts,1);

for Δe = 1 : Ne, % scan through each triangle element Δe Compute coordinates of the circumcenter of triangle Δe
for qi = 1 : 3,

% local node numbers qj and qk (e.g., if qi = 2 then qj = 3, qk = 1) qj = mod(qi,3) + 1;
qk = mod(qj,3) + 1;

% Get global node numbers ni, nj, nk ni = T(Δe,qi);
nj = T(Δe,qj);
nk = T(Δe,qk);

% Compute contributions to matrix elements from triangle Δe: Compute Lij and Wije for edge i-j
A(ni,nj) = A(ni,nj) + Wije/Lij; A(ni,ni) = A(ni,ni) - Wije/Lij;

Compute Lik and Wike for edge i-k A(ni,nk) = A(ni,nk) + Wike/Lik; A(ni,ni) = A(ni,ni) - Wike/Lik;

% ---- Contribution to RHS vector b from triangle Δe: Compute area Aicb of triangle Δicb
Compute area Aiac of triangle Δiac b(ni) = b(ni) - Re*Pz*(Aicb + Aiac);
end
end