BIEN/CENG 2310

MODELING FOR CHEMICAL AND BIOLOGICAL ENGINEERING

HONG KONG UNIVERSITY OF SCIENCE AND TECHNOLOGY, FALL 2024

HOMEWORK #3 (DUE OCT. 28, 2024)

1. In this problem, we will simulate the motion of a planet orbiting the sun:

As shown, we place the sun is located at the origin, and the aphelion (the point at which the planet is farthest away from the sun) is located on the positive -axis. At time , the planet is

located at the coordinate  and its distance from the sun is = √2 + 2. We can model the motion of the planet as a set of two second-order ODEs:

where  is the gravitational constant,  is the mass of the sun, and the combination  is equal to 2.94 × 104 AU3d2. (For this problem, we will use the time unit of days (d), and the length unit of astronomical units (AU), which is defined as the average distance from the Earth to the sun, about 149.6 × 106 km. ) We will choose the initial location of the planet to be the aphelion, namely, ( = 0) = 0, (= 0) = 0. We also know that, since the orbit is an

ellipse, at the aphelion | = 0. The speed of the planet at the aphelion is 

(a) Write a MATLAB program to solve the set of two second-order ODEs as an initial value problem. Allow the user to specify 0, the distance of the planet from the sun at the aphelion, and 0, the speed of the planet at the aphelion.

The program should stop when the planet returns to the aphelion, and output the period , the time it takes to complete one cycle. Your function definition should be:

                  function tau = solarIVP(R0, v0, showplot)

If showplot is set to true, provide a plot that shows the planet (a blue circle) moving around the sun (a red circle) as a movie. The speed at which the planet moves in the movie should be proportional to the speed it actually moves in orbit around the sun.

(b) Suppose we have a planet for which we can measure its distance from the sun at the aphelion, 0, and the period  of its orbit. Solve the boundary value problem to determine its speed at its aphelion 0, using the shooting method. Your function definition should be:

                       function v0 = solarBVP(R0, tau)

There is no need to produce any plot or movie for this part.

Hint: A good initial guess of 0 is √. You may call your function from Part (a). Some data to test your program (do NOT expect exact match):

    Planet

Mercury Earth Mars

0.46670 1.01673 1.6662

87.969 365.25 687.98

0.02269 0.01692 0.01271

                DELIVERABLES:

Submit your programs solarIVP.m and solarBVP.m. No need to provide any write-up or plot for this question.

2. To help cool down computer chips, heat sinks like the one shown below are often employed to carry away the heat generated more efficiently:

Consider one of the metal pins, represented in the following schematic diagram:

  Convection

where the temperature  is a function of both time and location (measured axially from the root of the pin),  is the thermal diffusivity that measures heat conduction in the metal,  is a parameter that measures heat convection from the metal pin to the surrounding air, and  is the temperature of the air around the pin.

(a) Suppose we are only interested in the steady-state temperature profile of the pin, i.e., when the computer chip has been running continuously for a while, and ejects a constant flux of heat to the pin. The PDE can then be simplified to a second-order ODE for :

       Hot computer chip at constant

Air at constant temperature