Hi,

I am trying to write a simple 2-D code (in FEM) to model a heat exchanger. The problem with my code is that it gives me some weird results.

The problem that I am trying to solve is that of two concentring tubes (the inner one is out of copper, and the outer one is plastic). To simplify my results, I am considering the 2-Dimensional configuration of the above problem, and hence I have the following configuration:

a slab with thickness do, length l (outer tube)

a slab with thickness dc, length l (the copper tube)

a slab with thickness di, length l (inner tube)

At the interfaces I impose continuity of heat flow and continuity of temperature. The free end of the outer tube is insulated, and I impose symmetric b.c at the free end of the inner tube. I consider the case of counterflow, where at the left end the temperature is To for the outer tube, and at the right end the temperature is Ti for the inner tube. For the rest ends, I assume zero flux.

When, Newton Raphson is converged (based on the residuals), the temperature in the copper has some really large values. I believe that this happens because of the several orders of magnitude difference in the values of the thermal coefficients. Which may lead to a stiff problem.

Is the above behavior normal or am I doing something wrong? If it is normal, how do I solve this problem?

Thank you

Vasilis,

Your problem definition is not fully clear to me: Do you consider only solids (i.e., no fluid flow in gaps between the pipes, and heat transfer is merely by conduction)? And what geometry, material and equations are modeled between the copper inner wall and the symmetry line/plane? Maybe a picture of the geometry, material zones and BCs would be helpful to clarify.

Independent of my doubts above, if you suspect your strange results are stemming from the large diference in the material properties, try to first solve the same problem with close (and even equal for the first shot) properties, and see if you have reasonable results.

I assume this is a steady-state pure conduction problem. What do you mean by "some really large values"?

Denote by Tf an intermediate temperature given by solving kc(Ti-Tf)=kp(Tf-To). kc, kp are the thermal conductivities of copper and plastic. Because of the huge difference in conductivities, Tf will be much closer to Ti than to To. Over most of the length of the two interfaces, the interface temperature should be roughly Tf. The temperature everywhere in the copper should lie between Ti and [roughly Tf], while the temperature everywhere in the plastic should lie between [roughly Tf] and To. Nowhere should the temperature lie outside the range [To, Ti]. If your temperature field violates these strictures, then you have some sort of numerical error on your hands. Your boundary conditions are okay, if they are correctly implemented. "zero flux" is the same thing as "insulated". The degree of departure from one-dimensionality will depend on the values of do, dc, di, but regardless of this, the temperature field should be as above.

Thank you for your replies!

Let me clarify my geometry:

At the bottom we have insulation, then there is a layer of thickness do where water flows, then there is a layer of thickness dc which is the copper, and finally there is a layer of thickness di where again water flows. At the top surface, I impose symmetric b.c.

The thermal conductivity of water is about 0.56 (S.I units) and that of copper is 400 (S.I. units). Because of that difference I am getting some really large values for the temperature. By the way, To=308K and Ti=288K. If I "increase" the thermal conductivity of water to 400, then I do not have the aforementioned problem. Nevertheless, the temperature in the top layer of water at the left end is smaller than Ti. Remember that I impose that the temperature is Ti at the right end, and that the temperature is To at the left end for the bottom layer of water. I believe that this awkward!

I considered the simplest case, i.e. zero flow in the two water layers, and the thermal conductivity of water is equal to that of copper! The Newton-Raphson converges in one iteration, but the temperature in the water of the top layer, and at the left end is larger than To, and it is ~315K. If I increase the number of elements in the x and in the y-direction, I find that the solution has not converged and that the temperature in the water of the top layer at the left end is approaching that of To.

What suprise me here is that although I am solving the simplest problem, the solution has not converged with the number of elements.

One more question: Does symmetric b.c. mean zero heat flux?

Well, that changes everything, doesn't it?!

I wish you had clearly stated your problem in the first post. Even now, not everything is clear. You still have not specified what "very large values" means. Anyway, your last test with zero flow and k_w = k_c shows that an unphysical overshoot of about 7K occurs on a coarser mesh. I suspect the difficulty comes from your boundary conditions of "zero flux" at the left end for the upper water layer and at the right end for the lower water layer. For pure conduction, to answer your question, the symmetry bc implied zero (conductive) heat flux. At the left upper and right lower, presumably cold and hot water streams are exiting your domain. What does "zero flux" mean here? Presumably you set the horizontal first derivative of temperature to zero. A better condition is to set the second derivative of temperature in the horizontal direction to be zero. Remember that the exiting water of both cold and hot streams are carrying away heat.

Hi Ananda,

I apologize for my first post. I had the correct picture in my mind, but I did not use the correct words to describe it. By the term "very large values" I mean something like 8000K.

At the moment I am trying to obtain an analytical solution of the problem that I described in my previous post. I should be able to obtain one, and use it to validate the numerical solution (which apparently is wrong).

I would like to clarify one thing: Symmetric b.c. means that the normal part of the conductive heat flux is zero, and not the tangential part of it.

Now, you did raise an important question that troubles me a lot. What b.c. should I use at the two ends? If I set the normal part of the conductive heat flux to zero, I assume that I either have symmetric b.c. or that that end is insulated. What is the physical meaning of zero second derivative?

In fluid dynamics we usually impose zero normal stress, and the physical meaning of this b.c. is that the fluid is in contact with air. Could I apply a similar b.c? Namely, could I assume that q=h(T-T_air), where q=-KdT/dn?

Now, there is another thing that bothers me. Since, I have fluid in contact with solid, at the solid-fluid interface I should impose q=h(T_w-T_s),where q=-KdT/dn, and T_w is an averaged value of the temperature in the water. But, the simulation does not converge with this b.c. and thus I am using the q_w=q_s b.c. instead. This is another issue that I need to solve!

Well, 8000K clearly indicates a numerical error. Heat is entering the domain and encountering great resistance in leaving the domain.

You have not explained about how you treat the motion of the water. Is it imposed? Or is it (pre)calculated? Is it inviscid uniform flow? Or is it viscous boundary layer flow?

Except for the case of zero motion of the water, you generally have both conductive and convective transport of energy by the water streams. Assuming that the flow is incompressible, the water motion (mass & momentum) can be solved for (or imposed) first. This is followed by solution of the energy equation in the whole domain as a separate step, treating the water motion as known. You can subtract out the mechanical energy balance from the total energy budget, so that you deal with only the internal energy equation.

The internal energy equation in the copper domain of course reduces to only conductive heat entering/exiting each face of each control volume in your mesh. In the fluid part of the domain, there are not only conductive fluxes at each control volume face, but also convective fluxes of thermal energy. I assume you have already accounted for all of this in your discretization. Unless the water is moving very very slowly, at the inflow and outflow of each fluid stream, the convective thermal fluxes will generally be more significant than the conductive ones. Are you upwinding the convective thermal fluxes? Upwinding is probably a good strategy. At the inflow of each stream, the effect of the convection is to bring in the temperature information from outside the domain, and therefore your strategy of imposing To and Ti at these boundaries is reasonable, provided that the energy equation (in the layer of cells next to the outflow boundary) retains the convective thermal fluxes to allow heat to be advected out of the domain by the water. At the outflow, a similar upwinding philosophy would seem to indicate that the interior temperature should be convected to the boundary, which would be consistent with the "zero normal conductive flux" condition (dT/dx=0) that you use. However, I think that a second derivative equals zero condition would be better. It basically says that the thermal profile of the stream is "fully developed" in the sense of fully developed boundary layer flow.

The attempt to use a convective heat transfer coefficient "h" at the fluid-solid interface is probably inappropriate if your water model includes the detailed simulation of the boundary layer. The heat transfer physics summarized by the convection coefficient is an encapsulation of the conduction that takes place in the boundary layer and the advection of energy by the fluid motion. The convection coefficient approach is therefore more suitable for inviscid flow.

One test you might try is to set the conductivity of copper to be very low, say 0.001. Then the two water streams should basically not interact with each other, and their outlet temperatures should be close to their respective inlet temperatures.

Because I am trying to solve the simplest case, I am considering a plug flow. Therefore, I am not solving the Navier-Stokes, and I am assuming that only the x-component of the velocity is not zero.

I am solving this set of eqns with FEM, and so I can not do the upwinding that you mentioned or better I do not know how to do it.

Finally, my model for the water does not simulate the thin boundary layer close to the solid surface. Because of that, I should use the convective heat transfer coefficient.

I see. Well, are you solving a set of differential equations, or just the energy conservation equation? In the energy conservation equation, are you properly retaining the advective energy fluxes at the outflow boundaries? If there is any problem with those, then heat will build up inside the domain leading to high temperatures large enough to drive the heat out of the domain at the cold stream INFLOW by conductive thermal flux working upstream against the advective thermal flux.

When you say "zero flux" at the outflow boundaries, do you mean "zero conductive flux" or "zero total (conductive + advective) flux"?

Notice that the natural way for the heat to exit the domain is at the outlets of the hot and cold streams. That test I mentioned with the k of copper set very low would help isolate the problem.

That is all I can think of. Good luck tracing the problem.

I did not quite understand the meaning of your first paragraph.

I set the conductive flux to zero. I will definitely try your test..