1-D Conduction+Convection: eriving a Finite Difference Scheme

[Saladsamurai]

1-D Conduction+Convection::Deriving a Finite Difference Scheme

***For those of you who would like to skip the derivation of the energy balance, skip ahead to the bold subheading below where the actual Finite Difference Scheme begins***

I am trying to derive a finite difference scheme for 1-Dimensional conduction and convection. I am using a similar approach to one used in my lecture notes for conduction only. It is a control volume approach to this.

I would upload a picture, but I think the 1D case is simple enough to picture. Picture a bar with constant cross-sectional area. We will work with a cross-sectional "slice" of the bar whose length is dx

EDIT: Picture added

So as said, there is conduction occurring (i.e. a temperature gradient is present) and there is convection going on. Here are the properties of the bar [with dimensions]:

A = constant cross-sectional area [Length²]
K = thermal conductivity [Power/(Length-Temp)]
c = specific heat capacity [Energy/(Mass-Temp)]
h = convective coefficient of heat transfer [Power/(Length²-Temp)]
\rho = density [Mass/Length³]The Energy balance is given by:

E_{in} + E_{generated} = \Delta U +E_{out}\qquad(1)

There is a conduction term on the left edge of the control volume at x and a conduction term at the right edge at x + dx

That is, the energy balance becomes:

q_xA\,dt+QA\,dx\,dt=\Delta U+q_{x+dx}A\,dt + q_hP\,dx\,dt \qquad (2)

where:

q_x is the heat conducted (heat flux) into the control volume (CV) at the surface edge x.
q_{x+dx} is the heat conducted (heat flux) out of the control volume at the surface edge x+dx.
Q is the internal heat source (or sink); heat generated (or removed) per unit volume.
q_h = h(T-T_{\infty}) is the convective heat transfer
P is the perimeter around the constant cross-sectional area A

Presumably they have neglected convective heat transfer (HT) at the end faces and are only considering HT via the body's surface area *though I am not sure why ~ Any thoughts on that?).

So here goes the derivation:

first we note that the conductive term at x is given by Fourier's law of conduction:

q_x = -K_{xx}\frac{dT}{dx} \qquad (3)

Using a two-term Taylor series, we can expand the q_{x+dx} term to be:

<br /> q_{x+dx} = -[K_{xx}\frac{dT}{dx}+\frac{d}{dx}(K_{xx}\frac{dT}{dx})\,dx] \qquad (4)<br />

plugging (3) and (4) into (2) yields:

-K_{xx}\frac{dT}{dx}A\,dt+QA\,dx\,dt=\Delta U-[K_{xx}\frac{dT}{dx}+\frac{d}{dx}(K_{xx}\frac{dT}{dx})\,dx]A\,dt + q_hP\,dx\,dt <br /> \qquad (5)

Canceling the -K_{xx}\frac{dT}{dx}A\,dt term of both sides and then dividing through by A*dx*dt and noting that \Delta U = c\rho A\,dT\,dx and assuming that K_xx is constant yields:

K\frac{\partial{}^2 T}{\partial{x}^2} + Q = \rho c \frac{\partial{T}}{\partial{t}} +\frac{hP}{A}(T - T_{\infty}) \qquad (6)

or for no source term we have:

\frac{\partial{}^2 T}{\partial{x}^2} =\frac{ \rho c}{K}\, \frac{\partial{T}}{\partial{t}}+\frac{hP}{AK}(T - T_{\infty}) \qquad (7)*****************************************************************************
*****************************************************************************
Now for the finite difference scheme:

If we have a one dimensional bar we want to find the temperature at a point 'p' that is located equidistant in the x-direction to points 'p+1' and 'p-1.' We can rewrite the derivative terms via Taylor expansions; that is,

\frac{\partial{}^2 T}{\partial{x}^2} \approx \frac{T(x+dx) - 2T(x) + T(x-dx)}{dx^2}<br /> \qquad (8)

and

\frac{\partial{T}}{\partial{t}} \approx \frac{T(t+dt)-T(t)}{dt}\qquad (9)

Hence, (7) becomes:

\frac{T(x+dx) - 2T(x) + T(x-dx)}{dx^2}} =\frac{ \rho c}{K}\, \frac{T(t+dt)-T(t)}{dt}+\frac{hP}{AK}(T - T_{\infty}) \qquad (10)

EQUATION (10) is where I am a little confused. I don't think that these should be "total derivatives" but instead should be "partials" but as engineers, we like to abuse notation. So I would like to fix this (any thoughts are helpful).

Looking at the right-hand-side of (10) there are 2 "T" terms I need to work out. Is the "(T+dt)-term" at location "x" ? I believe it is. And is the convection "T-term" at time "t" ? I also believe that is true.

If so, I am going to let T(x) at time "t" be represented by T_p^t similarly T(x,t+dt) would be represented by T_p^{t+dt} and so on.

With this notation, we can represent equation (10) as:

T_{p+1} - 2T_p + T_{p-1} =(\frac{ dx^2\rho c}{K})\, \frac{T^{t+dt}_p-T_p^t}{dt}+(\frac{dx^2\,hP}{AK})(T_p^t - T_{\infty}) \qquad (11)

Now we have purposefully left the time dependence off of the T-terms on the left hand side as this will determine whether the scheme is explicit or implicit --that is whether it is evaluated at 't' or 't+dt'

I am still unsure of the "T" on the right-hand-side in the brackets (T_p^t - T_{\infty}).
Should I have left the time dependence ambiguous on that as well? And define it when I define the scheme type? And is it indeed at point 'p' (pretty sure it is)?

 

[Mapes]

EQUATION (10) is where I am a little confused. I don't think that these should be "total derivatives" but instead should be "partials" but as engineers, we like to abuse notation. So I would like to fix this (any thoughts are helpful).

They're neither; they're finite times and distances \Delta t and \Delta x (thus, "finite-difference scheme ).

Everything else looks good. When you discretize a problem, you've got to make some decisions about which times will apply to which nodes in the approximation. These are arbitrary decisions; in the limit of small distances and time steps, all solutions will converge. If you already know about explicit and implicit solution methods, you're in great shape.

 

[Saladsamurai]

Hey there Mapes Ok, I think I see what you are getting at. I am trying to gear up for a FEA exam this week that will cover some FD methods too. Unfortunately, we did not learn FD from the text (because it was not covered), but instead from a couple of lectures. My professor presented the case with conduction only as follows:

<br /> T_{p+1} - 2T_p + T_{p-1} =(\frac{ dx^2\rho c}{K})\, \frac{T^{t+dt}_p-T_p^t}{dt}\qquad (12)<br />

He said that when we evaluate the LHS at time step t it is explicit and when we evaluate it at t+dt it is implicit.

So the left the LHS time dependence ambiguous. I am trying to extend the approach to the convection term, but I cannot figure out something:

Should my convection term have (T_p^t - T_{\infty}) in it or (T_p - T_{\infty}) in it and then base the time dependence on whether my scheme is explicit or implicit? I think it is the latter and here is why.

In the conduction-only equation (12) it is clear that the RHS is at the time-dependences given since the conduction term came from a Taylor expansion. However the convection term did not come from a Taylor expansion and thus I believe that its time dependence is ambiguous and should be scheme-dependent. Therefore, my equation (11) is not quite correct and instead should lose the stated time-dependence in the convection term.

Thoughts?

 

[Mapes]

The convection T_p is the same as the conduction T_p, so I'd recommend using the same time value for each.

 

[Saladsamurai]

Hey Mapes,

Unfortunately, the conduction term has both time steps in it since it is given by

(\frac{ dx^2\rho c}{K})\, \frac{T^{t+dt}_p-T_p^t}{dt}\qquad (13)

 

[Mapes]

That's not the conduction term, it's the transient or energy storage term. The k(d^2T/dx^2) term is the conduction term.

 

[Saladsamurai]

Yes, you are right, after I posted that I realized it did not make sense since conduction implies a Temp-gradient. So I have to work backwards for a minute here. OK! This is where the derivative thing was confusing me. Since my professor only gave us the finished equation, I am not sure where those 'T' terms ended up...give me a minute an I will figure it out and post back.

 

[Saladsamurai]

OK. I see it; that was easy. After dividing off 'K' we were simply left with the LHS of (12). So the time dependence for all conduction terms is ambiguous and thus the same should go for convection. So I should change (11) to

<br /> T_{p+1} - 2T_p + T_{p-1} =(\frac{ dx^2\rho c}{K})\, \frac{T^{t+dt}_p-T_p^t}{dt}+(\frac{dx^2\,hP}{AK})(T_p - T_{\infty}) \qquad (14)<br />

And perhaps I should move the T_p from (T_p - T_{\infty}) to the LHS so that all ambiguous terms are grouped together.

What do you think?

 

[Mapes]

I like it, but please make sure to write finite terms as \Delta x and \Delta t. Writing dx and dt is incompatible with the use of discrete spatial and temporal terms (e.g., T_{p+1}) elsewhere in the equation.

 

[Saladsamurai]

I like it, but please make sure to write finite terms as \Delta x and \Delta t. Writing dx and dt is incompatible with the use of discrete spatial and temporal terms (e.g., T_{p+1}) elsewhere in the equation.

Yeah. I was a little torn; I love my little differentials too much to let go of them, but the transition from x+dx to p+1 felt a little dirty.

Thanks for the tip! And for looking at this post. (I am thinking it might be a good one for the library that no one reads )

 

[Saladsamurai]

I'm not going to lie Mapes, I still don't like it. There are some things on the initial derivation that I am unsettled about; care to discuss? If so:

When you said

The convection T_p is the same as the conduction T_p, so I'd recommend using the same time value for each.

it made me go back and look at what we are doing. Looking at the differential element in the picture, we have 2 conduction fluxes. We were able, through some cancelation of terms, to relate both back to a single temperature gradient. But the convection term is kind of whacky. It relates the convective heat transfer via the whole element to the ambient temp and one single temperature inside of the element. But we just said that there must be a gradient else there would be no conduction.

Do you know what is going on here that I am missing?

~Casey

 

[Saladsamurai]

I'm not going to lie Mapes, I still don't like it. There are some things on the initial derivation that I am unsettled about; care to discuss? If so:

When you said

it made me go back and look at what we are doing. Looking at the differential element in the picture, we have 2 conduction fluxes. We were able, through some cancelation of terms, to relate both back to a single temperature gradient. But the convection term is kind of whacky. It relates the convective heat transfer via the whole element to the ambient temp and one single temperature inside of the element. But we just said that there must be a gradient else there would be no conduction.

Do you know what is going on here that I am missing?

~Casey

Anyone have any thoughts on this particular question?

 

[Mapes]

Looking at the differential element in the picture, we have 2 conduction fluxes. We were able, through some cancelation of terms, to relate both back to a single temperature gradient. But the convection term is kind of whacky. It relates the convective heat transfer via the whole element to the ambient temp and one single temperature inside of the element. But we just said that there must be a gradient else there would be no conduction.

There is one single temperature for each element (in the differential limit of zero element width, no gradient can exist within the element). The gradient only exists as a difference in temperature between any particular element and an elements on either side. Furthermore, these two gradients (the one between the center and left elements and the center and right elements) are not necessarily equal. If they are not equal, then the temperature of the center element will change with time (ignoring convection). I'm not sure if this answers your question?

EDIT: I should add, the conduction term is not a single gradient. It is the gradient of two gradients, the one on the left and the one on the right. That's why there's a denominator of (\Delta x)^2 rather than just \Delta x.

 

[Mapes]

Also, I just saw your question about convection on the ends of the bar. Although this isn't addressed by the differential element analysis and resulting governing equation, it is addressed by the boundary condition at either end of the bar, which would be, for example,

\dot{Q}=hA[T(x,t)|_{x=L}-T_\infty]

or, applying Fourier's Law,

-kA\left[\frac{\partial T(x,t)}{\partial x}\right]_{x=L}=hA[T(x,t)|_{x=L}-T_\infty\mathm{.}]

 

[Saladsamurai]

This is great info Mapes. Thanks for your time. I think this answers most of my questions for now, but I will probably revisit this thread with new thoughts.

 

[Mapes]

Any time! I'm glad you're learning the derivation of this type of equation. Many people finish a class in heat transfer and have no idea where the second derivative comes from, for example. The technique of Taylor series expansion is very useful in deriving approximate governing equations for complex physical systems. (For example, you can show that an energy minimum in any energy landscape, no matter how complicated, is parabolic for small perturbations, and can thus be modeled as a simple spring... but that's another topic!)

 

[Cyrus]

Great thread guys.

 

[Saladsamurai]

Any time! I'm glad you're learning the derivation of this type of equation. Many people finish a class in heat transfer and have no idea where the second derivative comes from, for example. The technique of Taylor series expansion is very useful in deriving approximate governing equations for complex physical systems. (For example, you can show that an energy minimum in any energy landscape, no matter how complicated, is parabolic for small perturbations, and can thus be modeled as a simple spring... but that's another topic!)

Me too! This is actually for a FEA course that I took out of sequence, so I take the "real" heat transfer course next semester. I feel like I will have a leg up on it, so to speak. Perhaps I will give this thread a revival when I go on summer break (if you can call 1 week a break :/) to talk about this whole energy landscape = simple spring thing. I think that I have an underlying idea of what you mean. The FEA equations for Heat Transfer and for Structural Mechanics are strikingly similar in that for both, there are Hooke's law types of equations and analogies. And these all arise as a result of applying the Theory of Minimum Potential Energy (if that's what it is called).

We'll talk.

~Casey

 

[Saladsamurai]

As a side note, things seem to get a lot trickier when the cross-sectional area is not constant. I thought it would be a simple matter of just throwing a factor of A back onto (6) but it does not appear to be. The conduction areas are now different and this makes things a little difficult. Also, the area associated with the internal energy term is up for grabs. I am currently doing a problem that looks like a trapezoid from the side, so the area varies linearly with x. Hence, I can get away with using the average area of the of region that bounds the point in question.

If I ever get through this solution, perhaps I will post it as a reference.