Microscale Heat Transfer Project Assignment

The analytical approach to determining bulk properties of an ideal gas based on molecular data. One key aspect of the analytical models is the assumption of an infinite number of molecules. Therefore, the analytical approach breaks down at the nanoscale, so computational approaches are used instead. In this project you will get a brief glimpse of some of the computational approaches used to modeling ideal gases computationally using the Direct Simulation Monte Carlo (DSMC) method.

1. One of the important aspects of kinetic theory is sampling from a probability distribution function f(x), where x is the property that is distributed among the molecules. In order to accurately model a system of molecules, random positions, velocities, and energies must be assigned such that the resultant distribution of properties follows that prescribed by the theory.

a. The easiest way to do this is to integrate the probability distribution function _-∞∫^∞f(x)dx = 1 since 100% of all sampled values of x lie between -∞ and ∞. This allows you to develop a relation for f(x). Then, the cumulative distribution function is defined as F(x) = _-∞∫^xf(x) dx. Setting F(x) equal to a random number R in the range 0 ≤ R ≤ 1 allows for accurate sampling using a random number generator.

i. Show that using this technique, a randomly assigned value of x for a uniform distribution of x between a and b (no values of x outside of the range a ≤ x ≤ b) is: x = a + R(b - a) [hint: f(x) is constant in the range a ≤ x ≤ b].

ii. How would you randomly assign a value of x for the distribution f(x) = cx, where all values of x lie between a and b?

iii. Use MATLAB to show plots of the distributions for 10,000 samples of x for both (i) and (ii) for a = 3, b = 12. Create a normalized histogram with 30 bins and compare the results to the theoretical curves of f(x)dx, where dx is the width of the bins.

b. There are situations where the method in part (a) doesn't work mathematically. Therefore, the acceptance-rejection method is used.

i. Show that the method in part (a) for the x-velocity component in an equilibrium gas f(v_x) = (β/π^½)exp(-β²v_x²) yields F(v_x) = ½(1 + erf(βv_x)), where erf is the Gaussian error function. This shows that you cannot solve directly for v_x by setting F(v_x) = R.

ii. The acceptance-rejection method is written as f^{^}(x) = f(x)/f_max, where f_max is the maximum value of the distribution function. Determine f^{^}(v_x) for the distribution in (i).

iii. In this technique, a value of v^{^}x is sampled using a uniform distribution between a = -3/β and b = 3/β. Write out the expression for sampling of v^{^}_x using the result from Part (a, i).

iv. The expression for v^{^}_x developed in (iii) replaces v_x in the expression for f^{^}(v_x) in (ii)-write out the final expression for f^{^}(v_x) in terms of R.

v. For a given randomly determined value of f^{^}(v_x) using (iv), the value is accepted if f^{^}(v_x) > R₂, where R₂ is a second randomly generated number. Run MATLAB for 10,000 points with this technique to show that it obtains the appropriate distribution by comparing to the Maxwell x-velocity distribution. Use β = 1 and 30 bins for collection in your histogram.

2. In this project we are only focusing on a 1-d velocity distribution.

a. Starting with Eq. 2.36, show that (v_x²)^- = k_BT/m.

b. If we assume that Eq. 2.62 is still valid for a 1-d distribution, then show that v_rms,x = 1/2β.

c. Combine parts (a) and (b) to show that β = √(m/4k_BT).

3. Our simplified 1-d DSMC simulation will test conduction between two plates at different temperatures. Consider a channel between two infinite parallel plates. The plate at x = 0 contains a temperature of T_c = 290 K, and the plate at x = L = 2 μm is at T_h = 310 K. The spherical monatomic molecules have a diameter of d = 0.2 nm and a molar mass of M = 30 kg/kmol.

a. Initialize a set of 10,000 molecules with randomly-assigned positions in 0 ≤ x ≤ L.

b. Apply the acceptance-rejection method to randomly assign x-velocity vectors for each molecule using a temperature of Ti = 300 K.

c. Divide the domain into N_cell = 25 cells in the x-direction. The cells are of equal thickness Δx. So, the minimum x-value for cell i is (i - 1)Δx, and the maximum x-value is iΔx.

d. Use a "for" loop to do the following:

i. Advance the positions of the molecules using x_new = x_old + v_xΔt, where Δt = 0.2 ns.

ii. Check to see if the molecules hit a wall by seeing if x_new < 0 or x_new > L. If this happens, then scale the x-velocity of the particle using v_new = -v_old√(T_wall/T_i), where T_wall equals either T_c(x = 0) or T_h(x = L). Determine the new position of the molecules at the end of the time step by examining the pre- and post-collision portions of the timestep.

iii. Randomly determine if collisions occur in the channel. For N_a molecules in the cell, then randomly select ∼N_a/2 collision pairs using the FindPairs function provided. The probability that a collision occurs in a given pair is calculated using

P^~ = (N_a/N_r)(gΔt/√2Λ_r)

where N_a is the number of molecules in the cell, N_r = ρ_rΔ_x, ρ_r = N/L, g = |v_x,1 = v_x,2|, and Λr = [√2πρ_r³d²]^-1. A collision occurs if a random number R < P^~. As a hint, use the floor command in MATLAB to help quickly identify which cell a molecule is in.

iv. If a collision occurs, then re-sample an x-velocity using the acceptance-rejection method to get a parameter v_coll. Use the temperature T_i in doing the sampling. The new x-velocities for molecules 1 and 2 are calculated using v_1,new = v_cm + v_coll, v_2,new = v_cm - v_coll , where v_cm = ½(v_1,old + v_2,old).

e. Run the simulation for 1000 steps, calculating the average the number of molecules ?N_a? in each bin. Plot this parameter versus the x-position of the center of the bins. Also, determine the percentage of molecule pairs that collide throughout the duration of the simulation. Physically explain any trend seen in the plot.

f. Repeat steps a-e for 1,000 molecules run for 10,000 steps. Comment on and explain any differences you see in the two sets of plots or collision percentage for the two cases.

g. Repeat step f for a molecular diameter of 2 nm. Comment on and explain any differences you see in the two sets of plots or collision percentage for the two cases.

Attachment:- Assignment Files.rar