讲解 CENG 2310、辅导 matlab设计编程
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:
𝑑2𝑥 = −𝐺𝑀 ( 𝑥 ) 𝑑𝑡2 𝑟3
𝑑2𝑦 = −𝐺𝑀 ( 𝑦 ) 𝑑𝑡2 𝑟3
where 𝐺 is the gravitational constant, 𝑀 is the mass of the sun, and the combination 𝐺𝑀 is equal to 2.94 × 10−4 AU3d−2. (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 𝑑𝑦| = 𝑣0. 𝑑𝑡 𝑡=0 𝑑𝑡 𝑡=0
(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 √𝐺𝑀/𝑅0. 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 𝑇 𝑎
Metal pin
𝑎
temperature 𝑇 𝑐
0
Conduction
𝐿 𝑥
Its temperature profile can be described by the following partial differential equation (PDE):
with the boundary conditions:
𝜕𝑇 = 𝛼 (𝜕2𝑇) − 𝛽(𝑇 − 𝑇 )
𝜕𝑡 𝜕𝑥2
𝑎
0=𝛼(𝑑2𝑇)−𝛽(𝑇−𝑇 )
𝑇(𝑥 = 0) = 𝑇 𝑐
𝑑𝑇| =0 𝑑𝑥 𝑥=𝐿
𝑑𝑥2
𝑎
where 𝑇 is the computer chip’s temperature, and 𝐿 is the length of the pin. (Here we are 𝑐
assuming that the “tip” of the metal pin is small compared to its length, so that the heat loss at the tip (in the +𝑥 direction) would be negligible.)
>
> > >
Solve this boundary value problem by the finite difference method, dividing the pin’s length into 𝑛 equal pieces. Your function definition should be:
function dTdx0 = heatSinkSteady(alpha, beta, Ta, Tc, L, n)
The program should plot the steady-state temperature profile 𝑇 vs. 𝑥, and return the value of 𝑑𝑇| . (This value is proportional to the maximum heat rate that can be carried
𝑑𝑥 𝑥=0
away by the heat sink while keeping the chip temperature constant.)
(b) Your model is helpful for designing a heat sink. Given that 𝛼 = 0.001 cm2/s, 𝛽 = 0.03 s−1,
𝑇 =300K,𝑇 =340K,whatvalueof𝐿(thelengthofthepin)wouldyouchoose?Explain 𝑎𝑐
your answer.
(c) Solve the PDE for the transient behavior of the heat sink (i.e. without assuming steady
state) using the method of lines. The initial temperature of the whole metal pin is 𝑇 . For 𝑎
the boundary conditions, this time, instead of fixing the computer chip temperature at 𝑇 , 𝑐
we will assume that the heat flux ejected from the computer chip is constant at the steady state value, i.e. the value of 𝑑𝑇| you get from running the program in Part (a). Stop the
𝑑𝑥 𝑥=0
program when it reaches steady state, and make two plots, a 3-D plot of 𝑇 vs. 𝑥 vs. 𝑡, and
a “contour plot” of temperature profiles at 10 different time points overlaid on the same plot. Your function definition should be:
function [] = heatSinkTransient(alpha, beta, Ta, Tc, L, n)
DELIVERABLES:
Submit your programs heatSinkSteady.m for Part (a) and heatSinkTransient.m for Part (c). For both, we will set up the discretization schemes and the boundary conditions in class, to help you get started.
Also submit the write-up for Part (b), which should come with a plot to justify your answer.
3. In this problem we will model the so-called “diffusion disc assay” for measuring the effectiveness of an antibiotic to stop bacterial growth. A small disc with antibiotic is placed in the center of the agar plate with bacterial culture, and over time, the antibiotic will diffuse outwards. If the antibiotic is effective, it will stop the bacteria from growing near the disc, resulting in an inhibition zone. An antibiotic’s effectiveness is defined by the concentration required to inhibit bacterial growth, called the minimum inhibitory concentration (MIC); the lower the MIC, the more effective the antibiotic is. In this assay, the size of the inhibition zone measured at a given time after applying the disc is used to calculate the MIC.
As shown, the agar plate is circular, and we place the origin at its center. The radius of the plate is 𝑅, and the radius of the antibiotic disc is 𝜀. At time 𝑡 = 0, we place the antibiotic disc, and the concentration of the antibiotic at 𝑟 ≤ 𝜀 is assumed to be constant at 𝐶𝑑𝑖𝑠𝑐 at all times. As the antibiotic diffuses outwards, the concentration of the antibiotic, 𝐶(𝑡, 𝑟), as a function of time 𝑡 and radial distance from the center, 𝑟, can be modeled by a PDE. At time 𝑡 = 𝑡𝑓, we measure the radius of the inhibition zone, 𝑅𝑧𝑜𝑛𝑒. The MIC is equal to 𝐶(𝑡 = 𝑡𝑓,𝑟 = 𝑅𝑧𝑜𝑛𝑒).
(a) By writing a balance equation for the antibiotic for the ring-shaped control volume on the next page, taking the limit of ∆𝑟 → 0, and applying Fick’s Law, show that the diffusion can be described by the following PDE:
𝜕𝐶 𝐷(𝜕2𝐶+1𝜕𝐶) 𝜕𝑡 = { 𝜕𝑟2 𝑟 𝜕𝑟
0
where 𝐶(𝑡, 𝑟) is the concentration of the antibiotic at time 𝑡 and radial distance 𝑟 from the
center, 𝐷 is the diffusivity of the antibiotic in agar.
State any assumption(s). Also, write down suitable initial conditions and boundary conditions. Note that for boundary conditions, it makes sense to impose “no flux” boundary conditions at 𝑟 = 𝑅 (the edge of the plate) and at 𝑟 = 0 (center of the plate).
𝑓𝑜𝑟 𝑟>𝜀 𝑓𝑜𝑟 𝑟 ≤ 𝜀
Hint:
lim (𝑥+∆𝑥)𝑓(𝑥+∆𝑥)−𝑥𝑓(𝑥)= 𝑑 (𝑥𝑓) ∆𝑥→0 ∆𝑥 𝑑𝑥
(b) To solve this PDE by the “Method of Lines” in MATLAB, we will divide the space domain 0 ≤ 𝑟 ≤ 𝑅 equally into 𝑛 pieces of width h, and call the concentrations at the boundary of adjacent pieces 𝐶1(𝑡), 𝐶2(𝑡), 𝐶3(𝑡), ... , 𝐶𝑛(𝑡), as shown below:
Important note: We will not simulate the point at exactly 𝑟 = 0, since it will result in a division by zero. Instead, we can assume that the concentration there is equal to 𝐶1(𝑡), in line with our “no flux” boundary condition. For our purpose, it will not matter, since that point will have constant concentration at 𝐶𝑑𝑖𝑠𝑐 anyway.
Let 𝐶𝑒 (𝑡) be the rightmost node within the disc, i.e., 𝑒h ≤ 𝜀 . Write down the ODEs for 𝐶1(𝑡), 𝐶2(𝑡), 𝐶3(𝑡), ... , 𝐶𝑛(𝑡), using finite difference approximations for 𝜕2𝐶⁄𝜕𝑟2 and 𝜕𝐶⁄𝜕𝑟.
(c) Complete the provided MATLAB program (antibioticDisc_template.m) to solve this PDE. Theprogramshouldacceptinputparametersof𝐷,𝑅,𝐶𝑑𝑖𝑠𝑐,𝜀,𝑡𝑓,and𝑅𝑧𝑜𝑛𝑒,plot𝐶(𝑡,𝑟) in a 3D surface plot, and return the MIC, i.e., an estimate as close to 𝐶(𝑡 = 𝑡𝑓, 𝑟 = 𝑅𝑧𝑜𝑛𝑒) as possible. Rename your program antibioticDisc.m and submit on Canvas.
For your testing, the following are some sample plots (with parameters specified in the titles). Note the asterisk marking the point 𝐶(𝑡 = 𝑡𝑓, 𝑟 = 𝑅𝑧𝑜𝑛𝑒).
DELIVERABLES:
Submit your type-written or scanned hand-written write-up for Parts (a) and (b).
Submit your program antibioticDisc.m for Part (c). To save you some trouble in plotting, you should start from the antibioticDisc_template.m provided to you, and only add your code where it is marked “% Add code here”.