These are draft notes from subtitles. Please help to improve them. Thank you!
When we want to model a wildfire, heat transport is obviously a vital phenomenon to consider. It's going to occupy us for the major part of this unit. To get started, let's look at this system. It's a tank that completely isolated from its environment. and that contains two different compartments one that is filled with water that's initially hot and another compartment that's filled with water that's initially cool. Now, I want you to describe the flow of heat from one compartment to the other and the temporal evolution of the temperature of the left part and the temporal evolution of the temperature of the right part. Temperature is just another expression for energy. Energy is being exchanged between those two compartments. To include that idea in our equations, we should be using the Kelvin scale for temperatures, not the Celsius scale and not the Fahrenheit scale. At least for an ideal gas, the temperature measured on the Kelvin scale is proportional to the energy content of that gas. To model the system, we could say that for instance that every hour 10% of the energy content of the left compartment flows to the right compartment and the same percentage per hour would flow from the right compartment to the left compartment.
What's the right model for this idea? Should the rate of change of the left temperature be 0.1/h times the difference between the right and left temperature? Should the rate of change of the right temperature be 0.1/h times the same difference should the rate of change of the right temperature be 0.1/h times the left temperature? Should the rate of change of the left temperature be -0.1/h * left temperature or should the rate of change of the right temperature be 0.1/h times the difference of the left and the right temperature? Check all that apply.
This differential equation would mean that T1 decays exponentially to zero. We're losing a certain fraction every hour. This does not make sense. With this differential equation, T2 would gain from T1. The red process would be in place, but we would never lose anything--this is not true. In the second differential equation, initially T2 is smaller than T1, meaning that this difference is negative. The rate of change of T2 would initially at least be negative which doesn't make sense T2 has to grow. The first one is correct. The percentage of 10%--0.1 per hour of T2 flowing into T1 and there is minus 10% per hour of T1 flowing out of T1 and in a similar fashion the lower one is true as 10% per hour of T1 flowing into T2 a plus sign here and there is minus 10% per hour of T2 flowing out of T2 minus. So actually what's of interest to us is only the difference of the temperatures. Footnote. So we could have been using degrees Celsius and degrees Fahrenheit all the way because we're forming this difference and are not interested in the total amount of energy contained in each compartment.
And now we come to a model that has way more compartments than the one before. A model for heat conduction in a wire. Let's draw imaginary lines after each millimeter of that wire and treat each millimeter as one compartment. Every single compartment is treated in the same way, so we can just look at one single of them to get an idea of what happening, let's take no. 8. Compartment no. 8 loses a specific percentage of its energy per time to compartments no. 7 and no. 9. To be specific, let's say it's going to lose 1%/ms to the left and 1%/ms to the right, but there will also an energy flow from compartment no. 7 to no. 8 with the same percentage of the energy of no. 7 and there will be an energy flow from compartment no. 9 to compartment no. 8. Again, the same percentage but of the energy content of compartment no. 9. So after a short amount (h) of time, the temperature of the compartment no. 8 will be its initial temperature plus that amount of time times 1%/ms to be gaining 1%/ms from compartment no. 7 and we are gaining 1%/ms from compartment no. 9 but we are losing 1%/ms to the left and to the right so we're losing thrice that percentage. One final thing to be doing to clean things up a little, let's get rid of this 1%/ms here. If we work with seconds instead, 1 second amounts to 1000 ms, so we need 1000%/s. 1000% is 10--so now, we have an equation for the temperature of the compartment no. 8 after one time step. Of course, this works similarly for all other compartments. We've just change the numbers. Its about the left neighbor, the right neighbor, and ourselves.
So now, we have an equation that can be implemented on the computer. Let's start by simulating the following situation. A wire is briefly heated by a candle so that some part of the wire is at the temperature of the flame and the rest of the wire is at room temperature is times zero We blow out the candle so that there's no more heating and then we watch heat conduction in action. How does this energy spread to the left and to the right. The code comes with some pretty fine constants for your convenience. It provides to erase one with the old temperatures at the beginning of the step and one for the new temperatures at the end of the step. And you have to compute the new temperature from the old temperatures for the compartment number i A note about the implementation in case you're wondering about this line here, we're exchanging the roles of these two arrays temperatures new, temperatures old after every step. Once the new temperatures have been computed, they become the old temperatures for the next step, which means that one of these arrays can be reused for the new, new temperatures if you will. This saves us from using tons of arrays. We're just using these two arrays with different roles.
The implementation is almost of a verbatim copy of the equation. We just introduced one little changed because the temperature of the old compartment is used twice. So we introduced an auxillary variable to store it. The diagram starts with the initial distribution of the temperature. You see the action of the flame and then the heat starts to spread out further and further. And, of course, as the total amount of energy has to stay constant, we're not heating anymore--the peak has to go down. We didn't treat the endpoints in our simulation. So after this final step, we should be careful about what's going to happen.
This strange expression appeared in our equation for heat conduction along the wire. Now, we look into its mathematical meaning. The scary name for this topic is central difference formulas for derivatives. This is about determining the derivative, that means the slope of the tangent of functions that are based on empirical data, measurements, and now your approach to determine the slope of this function at this point would be to look a little further to the right, determine another point of that function there, and then connect this point to our original point and you can see that this doesn't really work well but we want to determine the derivative of our function namely the slope of the tangent. The simple trick is to go a little to the left and a little to the right. Determine the points of the function, left and right and then connect these. To be clear, this works far better. Now we can write down some equations. That's Δx stand for the total distance from left to right. Δ always is denoting the difference. The estimate for the slope would be what we gain in height divided by how far we went to the right. The denominator is simple, Δx is how far we went to the right. The numerator is a little more complex. It's the difference of this value of the function minus that value of the function. That's what we are gaining in height, so it's our function on the right-hand side, f(x₀ + Δx/2) from x₀, we are moving by half of this distance to the right plus Δx/2 minus the value of the function on the left-hand side f(x₀ -Δx/2) that go into the left by half of the distance minus. Now you can get an idea about why this is called central difference. Obviously, it's a difference and the difference is centered at x₀. we're moving as much to the right as we're moving to the left. Now, we're going to apply the same formula to the second derivative. The rate of change of the rate of change, and think acceleration. We simply use this formula not for the original function, but for its derivative. You replace f by f' wherever you find an f. So this f becomes an f'. This f becomes an f' and this f becomes an f'. Prime, the second derivative. But now we can use our original formula for these derivatives. This is going to get a little long. When estimating this first derivative, we get an expression with denominator Δx and in the numerator, we got the following. We are moving to the right by Δx/2, so we're ending with f(x₀+Δx). We are moving one further Δx/2 to the right minus we are moving Δx/2 to the left which means this is f(x₀) and now for this derivative. We are moving by Δx/2 to the right minus Δx/2 + Δx/2, canceled. So this becomes f(x₀) minus now we have to move by Δx/2 to the left which becomes f(x₀ - Δx)/Δx. Let's now collect what we have found in total, we are dividing by (Δx)² and in the numerator, we've got plus the value one step to the left plus the value one step to the right minus twice the value at our original position. So what originally appeared in our equation, T7 + T9 -2(T8) is 1 mm² plus the second derivative of the distribution of the temperature.
Given this central-difference formula for the second derivative, try to figure out essential difference formula for the fourth derivative. The rate of change of the rate of change of the rate of change of the rate of change. Obviously, there's is a factor of 1/Δx⁴, but how much of these values of the functions do we need as we go one step to the left, two steps to the left, one step to the right, two steps to the right. Insert the correct numbers.
Think of this like a pattern. It's 1 times the value of one step to the left--1 times the value one step to the right-- minus twice the value at the original position. And now we have to insert this pattern into itself. If you look at this one, one step to the left, we have to go one further step to the left times 1, this 1 times -2 and this 1 times 1. So one further step to the left is this one, two steps to the left minus 2. This was the one, one step to the left--1 is our original one. Similarly for the right one, we go one further step to the right and so on and so on. One step to the left and minus 2 times the original one. For the central one, we have to multiply the pattern by -2, so we've got -2, 4, -2, which boils down to 1 times the value 2 steps to the left minus 4 times the value 1 step to the left 6 times the value at the center minus 4 times the value of 1 step to the right and 1 times the value 2 steps to the right.
We use this equation to describe the conduction of heat along the wire and over in a position to put some mathematical meaning to it. As time goes on, the value of the temperature changes by the length of the time step times this expression. So this expression has to be the rate of change of the temperature with respect to time. Ṫ₈ if you will at times 0. Of course, we're always dealing with estimates here even though I'm writing an equal sign. But now that we know about central-difference formulas, we have another interpretation of the term. It's 1 mm squared delta x squared times the second derivative of our temperature. Of course, I'm again cheating a little here. We only know the temperatures of the different compartments. We don't have any curve of which we could form the second derivative. Now, we're going to write this equation in a professional fashion. Let's do this right. Temperature depends both on space and time. It changes along the wire that is both position, and it changes with time. We see temporal evolution. What we need here is the rate of change with respect to time and the way to be writing this is this. It's called the partial derivative of the temperature with respect to time. It's not written with lower case letter d but with curly type of d. This derivative of the temperature with respect to time equals 10 over 1 second times 1 mm². And this one is the second partial derivative of temperature with respect to position written like this Note how this is being written. The curly d is squared not the complete expression in the numerator. Even though in the denominator, everything is being squared. They're going to discuss partial derivatives in the next segments so bare with me. This equation is called the heat equation. In this case, it describes the conduction of heat along the wire along one dimension, and it's a partial differential equation. PDE is the technical term. It contains partial derivatives. That's why the differential equations we have seen so far are ordinary differential equations, ODEs. PDEs, partial differential equations, are typical for problems in space or problems in space and time. The constant involved here is often called alpha the thermal diffusivity. The name already hints at this being a diffuse equation. The heat diffuses along the wire.
Now, we have a function that depends on two variables. The temperature depends on time as well as on position. The typical way to visualize this is to use a 3-dimensional drawing. For every time and every position, we go up by the value of the function, and so in the end, we get some kind of curve surfers that describes the function. The partial derivative with respect to time computes the rate of change as time increases, but we leave x constant. The slope of this tangent line is the partial derivative with respect to time. The partial derivative with respect to x keeps the time constant and only looks at what happens when we change x towards the slope of such a line. In the end, it boils down to computing regular derivatives and simply leaving all the variables fixed that are not mentioned. If you form the partial derivative with respect to T, we treat x as though it was a constant. If you form the partial derivative with respect to x, you treat T as though it was a constant.
Here comes a simple situation in which to apply the heat equation. Assume that the heat is distributed in such a way along the length of the wire, the temperature depends linearly on position--one end is cool and one end is warm. If you wait for a short time, what happens with the central part of that function. Will it move left, will it stay as it is, will it move right?
The correct answer is it will stay as it is. One way of seeing this is the following--the central compartment loses temperature to the left and gains temperature from the right if you will. The amount of gain and loses proportional to the difference in temperature, but to the left and to the right, the difference in temperature is the same, so we lose as much as we gain. The temperature of the central compartment has to stay constant. Of course, the same is true for the compartment to the left. We'll just look at one further compartment down the road and see that this compartment to the left gains as much as its going to lose. Another way of seeing this behavior is to look at the heat equation. The partial derivative of the temperature with respect to time equals thermal diffusivity times the second derivative of temperature. Now, with respect to x, this is a line of constant slope, which means we have to form the derivative of this derivative with respect to x, but when the first derivative is constant, the derivative of this derivative is zero--this line has a slope of zero, so the right hand side of the heat equation equals zero, which leads us to the statement that there is no change in time. In reality, of course, we can never have this type of linear gradient going on and on. The temperature must not shrink below 0 K, but of course in practice, there's some upper limit as well. So eventually, we will see some changes to that curve because energy accumulates downstairs and is depleted upstairs.
This is again the simulation of heat conduction along a wire. Now, thermal diffusivity and Δx² are spared out and only the result after 10 seconds is being shown--no intermediate results. The result after 10 seconds is going to be complex enough because we're going to modify the step size. Find the step size that barely works. We don't want to see temperatures below the ambient temperature and in particular, we don't want to see negative temperatures. Leave three decimal places of that step size.
The kind of method that we are using to solve the partial differential equation, the heat equation that is, is professionally called the finite-difference scheme, because we are forming differences of values that fit in finite distances not infinite decimal distances. Diagrammatically, what we are doing looks like this. They are using the current value of our current position and the current value to its left and the current value to its right to estimate the next value at our current position. There are other ways for doing this. This is the most simple, but obviously, it leads to instability. In the previous quiz, we saw that this scheme has a particular issue with this sort of function. A function that oscillates up and down with each step in space. So let's compute what happens with this Berry function. Assume that this is our initial distribution of temperature. Temperature at times 0 depending on x. Of course, physically speaking, negative temperatures won't make sense and particularly are not without the unit of measurement. But now, we're doing a little mathematics. Now, let's compute what's going to be the next temperature at x equals to 0. The result of the first step of the temperature at time h, h is the step size, and position 0 will be its former value times 0 and position 0 plus step size times the rate of change, and the rate of change is thermal diffusivity times the second derivative with respect to x, which we estimate 1 over Δx². The old temperature to the left minus Δx plus the old temperature to the right plus Δx minus 2 times the old temperature at 0. So what's that going to be? The old temperature at position 0 equals 1. Same here. The old temperature to the left equals -1 and that's to the right is -1 as well. This is 1-1-1-2, which makes -4 step size times thermal diffusivity divided by Δx².
Let's look at this result--after one time step, the temperature at position 0 equal 1 minus 4 times the time step times thermal diffusivity divided by the square of Δx. What we want temperature should diffuse to the left and to the right. Eventually, everything should ever reach out--in this case, to zero. So we want to pick the step size in such a way that our function eventually decays to zero as time tends to infinity. What choice of time step do we have, does the time step have to be less than 1/4 * Δx² divided by thermal diffusivity or twice that ratio or exactly that ratio or 1/2 that ratio.
It's the final option. The time step has to be smaller than half of Δx² divided by thermal diffusivity. Let me explain why--if the time step was exactly equal to 1/2 * (Δx)² divided by thermal diffusivity. This expression would equal to and our result would be -1. What would happen is that our initial temperature distribution would be reflected after the first time step where it was 1 before, it would be -1 after, and it's easy to compute that where it was -1 before, it would be 1 after. The next time step would undo this reflection and here the blue distribution again and so on and so on. We would see infinite oscillation. We would see oscillation but no decay. If age is smaller than 1/2 this ratio, this expression is smaller than 2. If age is smaller than 1/2 times this ratio, this expression is smaller than 2, we're subtracting less, so the result is larger than -1. Oh a new temperature, 1/2 to 1 times step at position 0 is not -1 but a little more than that, and of course, it's similar for the others. So you see that when it totally reflecting our function, we're losing a little, and with the next time step, we again are going to lose a little, so that's going to work out well. As soon as age is smaller than 1/2 this ratio, we're going to see decay, but to stop that oscillation during the decay, we indeed have to pick age smaller than 1/4 of that ratio, so that we're subtracting less than 1. This condition shows that our method is really poor. If you use a spatial subdivision that's finite by a factor of 2, we have to use a time step that's smaller by a factor of .4. Hopefully, this reminds you of Unit 3, the recovered implicit methods to deal with such instabilities. For partial differential equations, implicit methods maybe unavoidable.
It's about time to get even one dimension more, heat conduction in 2D. So, we're talking about temperature distribution that depends on time, on x and as a premiere on y. One way for visualizing this is a movie. Each frame of the movie depicts the temperature distribution in x and y at that precise instance of time, then we come to the next frame of the movie which shows us the next distribution in x and y and so on and so on and what we expect to see is diffusion. If there's only one hot spot at the beginning, it should gradually dilute over time. The heat equation in 2 dimensions does not contain any surprise. The partial derivative of the temperature with respect to time equals α times the second derivative of the temperature with respect to x plus the second derivative of temperature with respect to y. Note aside, this expression is often called the Laplace operator applied to T. Don't confuse this triangle with the delta. It's also often called the Laplacian. You can imagine the Laplace operator to measure something like the bend of a surface in 2 dimensions and there's a very simple finite-difference scheme for that differential equation as well. Now, we're looking at one point at the current time and the current values left and right, below and above to compute the next value at that central position. That next value is temperature after one time step at x and y and we can spell out our estimate. This is the value before times 0 x and y plus the time step times the diffusivity divided by the spatial step size squared and now we got lots of terms The values of our function to the left, to the right, below, and above all at times 0 to the left (x - Δx0₁), to the right (x + Δx), below (y - Δx), and above (y + Δx). Hopefully, it's not already confusing that I'm using Δx to also denote the increment of y. Now that we've covered the outer point, the inner point remains. In one dimension, we are to subtract twice the value of the inner point. In two dimensions, we have to subtract four times its value. It's easy to see that this has to be a factor of 4. Imagine that it's equal to 100 Kelvin everywhere then we get 100 + 100, (+), (+) makes 400. We subtract 4 times 100 so we get 0. The new temperature will be old temperature as it has to be. There's nothing that could be smooth out if the temperature is already constant in space. This factor out 4 also hints at what's going to happen with instability. We have to be even more picky about the step size. It can only be half as large as in the one dimensional case. And the other artifacts that we are going to see close to instability are checkered board patterns.
Now, it's time to implement the diffusion equation in 2D. Of course, the arrays for the old and the new temperatures are now 2-dimensional arrays, and we initialize the arrays by putting in a right angle of high temperature. You may wonder about the order of y and x here. Shouldn't x be the first one and y be the second one. We use this order with the y appearing first and the x appearing last because we stick to the idea of the arrays being matrices. With the matrix, the first index refers to the number of the row and the second index refers to the number of the column. Changing the row mean to go up or down, which is like changing y, and changing the column means to go left or right, which works like changing x.
Once you've got the order on y and x right, the rest should be pretty straightforward. We're going one step left, we're going one step right, we're going one step down, we're going one step up, and subject four times the original value. And the result is that our initial square has turned into a soft blob. In case you're wondering, we're using 50 steps from the left to the right and from the bottom to the top. How could we get such a smooth picture with just 50 * 50 values. The trick is interpolation. So let's switch off interpolation to see how our raw data really look like. Now, you can see those 50 steps in x direction and y direction. Interpolation turns this pixelated image into one that looks really smooth.
Now, let's go shopping. What do we need to build a simple model of a wildfire. The most important simplification that we are introducing is to treat the forest as though it was a football field. We neglect height, we've reduced everything to two dimensions. In two dimensions, we can already describe how heat diffuses, but in addition to that, there's also heat loss, there's wind, and first and foremost, there's combustion. Now, we are going to look into these three starting with heat loss.
Which of these differential equations should we be using to model heat loss. We're introducing two constants here, the ambient temperature of 300 K and time constant for the heat loss of say 2 minutes. Should the rate of change of our temperature be that temperature minus the ambient temperature divided by the time constant or just our temperature divided by the time constant, or the difference the other way around divided by the time constant, or just the ambient temperature divided by the time constant.
Three of these solutions won't work already for the reason that they would lead to an increase in temperature even though the forest is warmer than the ambient temperature. With this solution, the rate of change of our temperature would always be positive. That doesn't look like heat loss. The same happens here. If the temperature of our forest is larger than the ambient temperature, we again get a positive rate of change here. This does not make sense. This is the right one. Think about two compartments. We have one compartment with temperature capital T and another compartment with a temperature of 300 K. The difference between these two is what controls the decay of the temperature, and of course, this difference is negative as it needs to be for our temperature to decrease.
Modeling the effect of wind seems to be pretty straightforward. It looks as though we simply have to push the fire through the forest with a certain velocity. If the left curve is the temperature distribution at times 0 after 1 time step, h will be shifted to the right curve, and the amount of shift amounts to the time step times the velocity. As we are storing these values in an array, the wind will simply shift their position in the array. The tricky thing, however, is that the time step times the velocity is really small. We would have to shift these entries by a fraction of a step. How would you do that. Here comes a different approach. I want to get the new value, the value on the right curve at the position X₀. I get that new value by looking it up on the original curve step size times velocity to the left from here, and now the trick is to use the derivative of this function. Let's write down what this new value is. The new value after 1 step time equals the time step h. At this position, x₀ is approximately the old value T(0, x₀) and now we have to take care of the difference between these two, and if we estimate it using the derivative, so all our correction the length of that green line will be step size times velocity this length times the derivative--this slope of our tangent. The derivative is the partial derivative of temperature with respect to x at times 0 and position x₀, and we have to be a bit careful about the sine. If the derivative is negative such as in this case, we have to go up, which means we need a minus sign in front. So that's the trick, we use the spatial derivative to accomplish this shift. As usual, however, as with all these methods, we have to be careful about not taking too large steps. This value with step size in time times the velocity has to be reasonably small. In particular, it should be pretty much smaller than our grit size Δx. The way to compute this derivative is, of course, to use a central-difference formula.
Now, implement the shift of the curve. There's no diffusion, no heat loss, no combustion. It's just about shifting the curve. The velocity is already specified and we have different initial conditions. Now, it's just the sine wave. We'll experiment with that in a second. You can already have a look at what's happening at the boundaries. They are pretty ugly.
No surprise in the implementation. We're subtracting step size times velocity times the estimate of the derivative and that estimate stems from our central-difference formula. In the residing diagram, you see that overall you achieved the right effect, the sine curve is being shifted, but we definitely need to take care of the boundary-- at least when our signal reaches the boundary. And let's try an experiment, I'm inserting the number 30 here, meaning that I want to start off with 30 periods of the sine--not the single period of the sine. And to make the image less clutter, let's use the number of 500 down here. This is the result that we get--we can see the artifacts at the boundaries, which we already know, which we already know, but there's a second phenomenon. We can see that the amplitude is slightly growing. We're growing from the blue curve to the green curve. So this looks unstable. The higher the frequency that we fit in, the more unstable this becomes. This is why I start with the sine wave and not with some rectangular function. The good thing is if we combine this way of shifting the function with diffusion, diffusion washes out the high frequencies and saves us from instability
It's time for the final ingredient, combustion. Here, we're just looking at one-singled cell at a time. We're treating each cell on its own. We're not looking at the neighbors unlike as we did with diffusion and with winds. For this singled cell that we're currently looking at, we have to take two different quantities into account. The first is the temperature which we already know, it's measured in Kelvin. The second is the amount of woods that can be burned. Actually that should be something like a density of wood to be measured in kg/m². As the combustion goes on in the cell, this density is reduced possibly to 0 and in addition we can use this density to model a variation in vegetation. In a thick forest, the value of this density goes up and at a location at which there's a clearing, the value of this density goes down. Now, we have to distinguish between two different cases. First, the temperature is below the ignition temperature which is something between 500 and 600 K. Second, the temperature is above the ignition temperature. Let me put an equal sign in here to really cover all possible cases even though in reality of this temperature, they'll never be exactly equal to that ignition temperature. When we are below that point of ignition, the amount of wood does not change. It's not burning. No heat is being produced because there is no combustion. Later on of course, the temperature is going to change through diffusion, through wind, and through heat loss but at the moment, I'm just looking at combustion. In the second case, we do have combustion. Wood is burning so the rate of change of wood with respect to time has to be negative. That's called the burn rate and the simple model for the burn rate would be to say that it's proportional to the current amount of wood. The more wood we have, the more we burn and we introduce a time constant for the burning that's called tB The larger this time constant is, the smaller the burn rate is, the less we burn. This time constant describes so to speak the time it takes to burn a specific amount of wood. Now we have to think about temperature. In the second case, wood is burning so the temperature should be increasing and that should be proportional to that burn rate. We need some proportionality constant in here versus something like the heating value converting from an amount of wood to an increase in temperature. This is our simple model of combustion. Now let's add the effect of heat loss in all equations about the rate of change of the temperature. We need to include heat loss which is + T amb-T/tHL for that heat loss and the same when there's combustion going on plus Tamb-T/tHL
As always, it's a good idea to look into equilibria, to understand what's happening. Most probably the time constant that governs the burning of wood is much higher than the time constant that's responsible for the heat loss. Heat loss will be quicker than burning. We can safely assume that the density of the wood is almost constant on the time scale on which the temperature develops, and now we can draw this nice diagram that shows how the rate of change of the temperature depends on temperature. If we are below ignition, we only have heat loss. That's the linear function with negative slope. And when the temperature equals the ambient temperature, the rate of change will be zero. So that's this part of the curve. A linear function with negative slope and that function is zero at the ambient temperature, and the second regime, when there is combustion, we have to add this constant. And as we assume that the density of the wood stays almost constant, the burn rate is almost constant. What we get is simply a shift upwards. The same linear function as before but shifted upward like this. Now, here's a question for you--how many stable equilibria do you see here. Enter that number.
It's relatively easy to get some reasonable values for the time constants, but it's pretty hard to estimate the heating value. There's just something that we can do with this curve. The reasonable value for this upper equilibrium is 1200 K. Just from the color of the flames, the wood should be burning at this temperature. Now, let's say that the ambient temperature is 300 K, the time constant for the heat loss is 2 minutes, the time constant for burning is half an hour, and we've got a wood density of 100 kg/m². From this data, compute the value of h--the heating value and enter that as an integer number of K/kg/m².
The trick is to add these two lines to the diagram. We know the slope of this line, which is -1 over the time constant for the heat loss. We know the length of this line, which is the jump, and the rate of change of the temperature. The heating value times the burning rate. This term makes the difference between burning and not burning, and now we can use the slope to express the ratio of this leg to this leg. 1 over the time constant for the heat loss equals the ratio of the vertical length to the horizontal length. The vertical length, the green line, is h times the burn rate and the burn rate is the wood's density over the time constant of burning, and the horizontal length is 1,200 K minus the ambient temperature. You solve this equation for h, plug in these numbers, and then you get 135 K/kg/m².
Now, we have all the ingredients to that model of a wildfire. The heat equation in 2D to describe the diffusion of heat. We have a model for heat loss. We can describe the effect of wind at least in one dimension, but that would be easy to generalize to two dimensions, and we have a simple model for combustion. In the homework, you're going to put all of these together to model a wildfire. This boils down to adding the different contributions to a rate of change of temperature. This is what the result could look like. In the region where the fire has started, the wood is depleted so they can't be any combustion going on. The fire has advanced to the right, a little to the left, and to the top.
In this unit, we look into spatially extended phenomena governed by partial differential equations rather than the ordinary differential equations that we have seen so far. Equilibria, again, played a vital role and we saw that solving partial differential equations on the computer leads to severe numerical issues. Next time, I'm going to give an overview on more advanced topics including some hands-on experiments.
In this problem, we're going to deal with heat conduction along the wire. Just like you saw in the unit, we're going to split the wire that we had into discrete chunks. In fact, we're going to split it into 100 different chunks and went from 0 to 99. We can then measure the temperature of each chunk and labeled that with a proper number. So for example the temperature in the last chunk is going to be labeled to 99. As you already heard from Jörn, the temperature in any given segment--let's take the 8th segment that's simple--depends on the temperature in the surrounding two segments. In those segments that are touching one another are going to feed energy into each other. So this segment and that segment can both be energy into the 8th. In the same way, the 8th segment is going to feed energy into both the 7th and the 9th segments. To express the temperature of the 8th segment after time step done, we can use the finite difference scheme to come up with this equation right here. This equation tells us that the temperature after the time step, right here, is equal to the temperature before the time step plus this term, which we can see depends on the initial temperatures of the 7th, 9th and 8th segments. This is obviously an explicit method for calculating T₈(h), but just for a change of pace since this is a pretty unstable solution, we're going to try an implicit method instead. Remember that implicit means that you're going to have this term over on this side as well. So we're just going to switch all these 0s to h's. We're going to use this implicit method, of course, converting it to code in this problem. One problem that we have to deal with is how to end with segments in other end of our wire. Both of these have an adjacent segment to one direction but not to the other. Since the segments start at 0 and count to 99, our problem was explicitly is that we don't have a segment numbered -1 nor a segment number 100. To deal with this, we're just going to set the value of T₋₁ to equal T₀ and similarly for T₁₀₀ set in out to T₉₉. This is going to be true for all times t. Your main job in this problem is going to be to define the coefficient for every equation in the set of equations right here. Each equation involves an expression, which is just the sum of all the different temperatures of every segment along with some coefficients stuck in front of certain one, and all these equals the initial temperature at that certain segment. There seem to be the set of a 100 different equations. If we're going to consider this system of equations as a whole, we can actually think of it more conveniently as dealing with matrices. This website right here numpy.matalg.solve deals with equations involving matrices. We have a matrix as the coefficients which are all these pink question marks multiplied by a single column factor, which contains all the temperatures after time step each. These multiplied together equals this vector, which contains all the initial values of the temperatures. As we already know the temperatures at times 0 and you have created a matrix of coefficients, we can use this method numpy.matalg.solve and input the main information and solve the vector that we want--the one containing the temperatures after time step h. Taking a look at our supplied code, we can see that we've created an initial temperature distribution for you. We've also created a parameter for you called c, which I think if you look back at the finite difference formula, you'll see if's going to be very useful for you. As I told you earlier, your first task is going to fill in the array named coefficients with the appropriate values. Remember, a lot of the coefficients will be 0. Another hint that I want to give you is that using for loops could be a great way to fill in your coefficients to all the different slots in the matrix. Once coefficients are filled in, come back up to the step size and set it to 0.5 and then submit. In the end, you should have a curve that looks nice and smooth like this despite the pretty large step size that we're going to be using.
Now let's talk about the solution to our problem with heat along the wire. First things first, I've rewritten the implicit equation that I gave you in the intro video. The only change here that you'll see is how I've replaced the coefficient in front of the parenthesis with C, the parameter that we created for you in the code. Now, if we rearranged this equation a little bit, we can get it in the form that's super convenient for using numpy.matalg.solve. This is because we already know T₈(0) and now we know what the coefficients are. So we are to solve for all of the temperatures after times that of h. This shows us then that two of the elements in the matrix coefficient are going to be -C and 1+2C. The one exception to this is going to occur at the boundaries. I told you these were going to be a little bit tricky, so let's talk about how to deal with them with right now. We take this equation right here but instead plug in the correct values for T₀. We see something interesting. There's a dependence on T₋₁. Now we defined our boundary with additions earlier, that T₋₁ is going to equal T₀ for any time. Once we plugged this n, we're actually going to get a pretty simple equation something that only depends on T₀ and T₁. To help you better visualize how are we going to fill in the array using loops, I've created this diagram right here which color codes the different possible coefficients that we could have in any slot in the matrix. As you just saw, a bunch of the terms in the matrix are going to be 0, so I filled in big 0s for those portions of the matrix. We know that the temperature of the second that we're looking at, is going to have a coefficient of 1+2C in front of it corresponding to this blue spirals. We know that both adjacent segments who have coefficients of -C which are these pink spirals over here. The two exceptions of course occur at the boundaries, for we have one pink spiral and one green spiral. For one coefficient -C and then one coefficient of 1+C. Now if I felt like working really hard, I guess I could have drawn this diagram with 100 rows and 100 columns making up the actual size of the matrix that we're looking at. But I was feeling lazy so I scaled it down a bit. However, if we think about it, we know that the number of blue spirals in that imaginary expanded diagram will be equal to 98 because there's going to be one blue spiral in every row with the exception of the first and last. Similarly, we would see 198 pink spirals, two for every row, except only one in the top and bottom rows. In the picture I just showed you, the blue spirals corresponded to the coefficient 1+2C and we said that we would need a loop using 98 iterations to express that. That's exactly what we have here in this top for loop. You can see that the next for loop contains one more iteration, since it starts at the index 1 lower. As expected, this corresponds to the pink spirals in the diagram with the coefficient equals -C. Here we have two exceptions for the two different boundaries and you remember that these were the green spirals in the diagram. We can check and see that we've set the step size equal to 0.5 seconds as we said we want to do. We can run the program and see what we get. Sure enough, we end up with the same solution that I showed you earlier. Since in our initial situation, we had a flame placed under the center of our wire, it makes sense that after some time has passed, we still see a peak in temperature toward the center of the wire, but the distribution has smooth out quite a bit. This has right in line that what we expect from the heat equation. Right now, the end time is set to 10 seconds which means that the final graph that we're looking at is the distribution of the temperature at this time. Let's see what happens if we set it to something smaller showing what's happening earlier on in the temperature change. So here's our distribution after just 4 seconds have passed. Consider we have a nice smooth curve, a lot like what we saw earlier except that this time the peak is at a much higher level than it was earlier. When end time was at the 10 seconds, the peak was just under 500 Kelvin and now we can see that it's resting right above 600. The distribution is also quite a bit narrower showing that we're farther away from the equilibrium position than we were with the larger step.
So now we come to the last problem of unit 6 and actually the last problem for entire course aside from the final exam at least. Now you have a chance to put together all of the informations that we gave you in unit 6 to create a 2-D wild fire. We're going to be looking at a big grid of forest and seeing how diffusion, heat loss, wind, and combustion come together to affect the time derivative of the temperature. As I talked about before, these four factors can simply be added together to show how T changes with time. Let's look at the code for a second to figure out what are we going to ask you to do. As you can see, they've given you a number of different constants that are going to help you out a lot. One thing that's kind of a small detail that maybe interesting to know is that the velocity due to the wind, which we can see right here, is actually very low. It's much lower than the velocity of the wind itself since the fire is only going to follow the wind at a fraction of the wind's real speed. As I said earlier, your first task is going to be to take into account heat diffusion, heat loss, wind, and combustion in the 2-D heat equation right here. Fortunately, we're not going to ask you to use the implicit method that you used in homework 1. You can just use the explicit method you used before so basically the forward Euler method. And important thing that we need to deal within this problem is especially the varying conditions. We must be careful to not introduce any hard edges anywhere, or they will probably lead to artifact. First in the solver to the huge values in the derivative and second in terms of the model itself. We want to make sure that we don't see any staircases or pixelated contours or really anything that does not look like reality. To help deal with this, we have already included a smooth calcium bell shaped for the initial distribution of the temperature. The second thing that you need to do, however, is to create an initial distribution for the wood density. Now, I'll give you a little bit more information about that right now. Remember the grid that we're creating using the code is basically making a map of the terrain that we're looking at. This graph right here is going to give you an idea of what we are looking for in terms of wood density distribution. Now here we have a bunch of lines that are parallel to one another but slightly displaced. X and y are actually just the physical coordinates that we are looking at. But the different colors of the different lines correspond to changes in the wood density at each location. So every point along the blue line is going to have the same density as every other point on that line. Every point on this purple line has the same density and so on and so forth. In the code, we created 2 constants for you, wood 1 and wood 2 which have underscores in the code. That would be very clear what they are referring to. In this mark out, the highest and the lowest possible values of wood density that we have in our forest. Now, all of these lines are parallel to one another, and we've given you the slope and the code. That means that the only factor differentiating from one another is their y intercepts. If you make a graph of wood density versus y intercept to make something looks like this. You can see that we are only letting wood density vary between wood 1 and wood 2, right here. Another hint that might be helpful when you're doing this problem is that you can solve for the y intercept of given line by saying that the y intercept equals some y value minus the slope times the x value of that point. So it's kind of a tricky problem but remember it's your very last one before the final. If you have any trouble, come to forums. Good luck.
Your first task was to fill in the rest of the code of this for loop down here. So the first thing that we changed is the wood density. The new wood density is going to be equal to the old wood density minus the time step times the burn rate. Since remember the burn rate is really just the time derivative of the wood density. Now for the change in temperature. Of course, we start with the initial temperature. And then we add to it h times a bunch of stuff. Now all of this code is just taking into account those four factors that we were saying influenced the change in temperature. So the first thing we have right here you can see is a 2D heat equation. And this is actually implemented using the finite difference scheme. After this, we also take into account heat loss and then the wind and then also the combustion. Now going up to the top of our wildfire function, let's talk about this initial distribution of wood density. This first line right here is just a translation of that mathematical expression I had written out for you. The y intercept of the point that we're considering is just equal to some y value minus the slope that we're given times x and actually this xy is the xy that we're talking about at this moment. That means that the wood density then is equal to wood value 1. Well, a bunch of typos. Well that's a problem. Let's fix that right now. See everyone makes mistakes. This should be wood 2. This should be wood 1. Okay now that all of that is fixed we have the wood density at our point xy is equal to wood 1 the maximum possible wood density plus all of this. This is the difference between wood 2 and wood 1. It's that range in which w is allowed to vary across the terrain divided by the difference between intercepts. So you remember on the graph that I showed you earlier. This graph right here actually. That is the slope of this line since wood density here is on the y axis and the intercept is actually on the x axis. If we take the difference in wood density divided by the difference in intercept we get the slope of this pink line right here. So now you can see that we have the slope of that pink line times the intercept of the line up to the point we're considering is on minus intercept 1. So this is really just using point slope form to calculate w. Remember, we only want w to range between wood 1 and wood 2. So this line restricts it to that range. Now let's look at the plot we get. So finally a graph of our forest. In both images, we're looking at the xy square of terrain that we've been talking about. Here we're graphing the temperature shown by these different colors, and here we're graphing the density of wood. Now unfortunately this is only a 50 x 50 pixel image. If we wanted to see what's happening with much better resolution, we should use 200 x 200. Unfortunately my computer is not cooperating so we're going to have to settle for this less perfect image. Even though this isn't really a perfect picture, you can still tell that the wood has been depleted in the center. The density of wood is very low right here compared to the rest of the area where the fire hasn't spread yet. You can also see that as we wanted the density of the wood does vary along lines like this. In the 200 x 200 image, you would be able to see that the wood had been very depleted in the center and that the fire had marched out the rim of this area. We can kind of tell that there is sort of a red circle around here. However, we'd be able to better that the fire is moving slower, and it's less hot in the area where there's less wood. The spreading the fire is way to slow on this image and one reason for that is that the heat is averaged over larger pixels so actually not lowering the temperature reached the temperature of ignition in a lot of places. I hope you found talking about wildfires as exciting as I have. Awesome job on this problem and throughout the entire course. I hope you're excited for the final exam. I think it'll give you a great look at a ton of different examples not just from the material we've learned about but also from some new kinds of applications.