Dear numerical expert,

I have a false numerical thermal problem.

There is a duct (2D) with length of long (20 times of height of the duct) with heated on top and bottom wall and with a inlet with cold fluid inflow (for example. temperature is zero and hot walls are 1.0).

When I applied a upwind difference scheme to the energy equation, the temperature in the duct seems reasonable but its accuracy is low as we know that at last it has first order acuracy.

To get a high order accuracy answers, I used the central diffence scheme in discretizing the convection terms. I got a converged answer but the temperature near the inlet showed minus solutions. Of course, it is not correct but it has the second accuracy answer!

By testing the blending factor between the upwind and the central diffence, the reasonable soltions rest in the upwind side. I have no available imagination to check if there are such a kind of solutions for the velocity component but it could be done for the energy transportation.

I think this kind of problem is not new but it really needs to solve. Central-Upwind hybrid sheme, Exponential scheme, QUICK, LOAD, CONDIF schemes exist but which one is the best, with high order spatial accuracy and with stable characteristics?

Is it correct by setting some limiters in the solving the eneergy eqs. but the central difference scheme is still adopted. To avoid the divergence, the upwind scheme is applied to the coefficient and the source terms due to the central diffence are modified.

I am expecting any advices.

Yun Mo

Hi.

Some gurus told me that Leonard's is the best. [Leonard B.P. A stable and accurate convective modeling procedure based on quadratic upstream interpolation. Comp. Meth. Appl. Mech. Eng., 1979, v.19, No.1, 59-98] It works well indeed, I tried.

However, you should be aware of certain complication in your particular problem. Your boundary conditions are in a sense incompatible. Namely, at the corners of the computational domain you have a singularity of the solution, since temperature is zero at the inlet and unity at the wall. The solution exists, but the derivatives tend to infinity as the corner point is approached. Therefore, the finite difference scheme does not approximate the solution near this point. To avoid oscillations you have to use a so-called monotone scheme, almost invariably of the first order.

There are two ways to deal with it, I tried both, and both work well. First, you can use a combined scheme. For example, you approximate dT/dx with

(1-a(i))*(T(i+1)-T(i-1))/(x(i+1)-x(i-1))+a(i)*(T(i-1)-T(i))/(x(i-1)-x(i))

With a(i)=0 you have second order central difference, and with a(i)=1 first order upstream difference. Now relate a(i) to the grid step so that for large steps a(i) is close to 0 and for small steps a(i) is close to 1. Last, use non-uniform grid x(i) such that at the inlet x(2)-x(1) were proportional to 1/N**2, N being the number of grid points, but far from it take x(i+1)-x(i) of order 1/N, and smoothly variing in between. You can construct the formula yourself. In that case you will have 1/N**2, that is second order, scheme and upstream differencing near the inlet. It works quite well.

The second approach is to exclude the vicinity of the corner from the computational domain and use asymptotics of the exact solution in that vicinity instead. This second approach is more complicated and requires certain asymptotic skills, really. And the order of accuracy is usually higher than first but less than second.

However, are you sure that you really need such boundary conditions?

Hope this may help. Yours Sergei

Dear Sergei,

I very appreciate your advise. I might not have given appreciative explanation. The finite volume method is adopted here. The model and the part of solution of steady temperature field in the upper region are shown. It can be seen that there is a minus temperature region near the inlet. The solution is converged. Though the mius value is not so large compared with the maximum temperature (=1.0) but it may not be allowed.

The values at 4 corners are not used in FV method. Near the inlet, both convecton and diffusion should be considered.

By switching the upwind scheme it seems correct but in a lower accuracy indeed. The following solution seems not completely correct but really in a second accuracy!

Other method such as QUICK method is somewhat complex in the diffencing and I am not used at present. I like simple but higher accuracy, easy use and check method. Too complex scheme may lead to higher possibility of wrong results.

Thank you very much, Sergei.

-------- T = 1.0 ------------------ upper wall

.....................................

----> T=0.0

----> U=0.1........................---> outlet

----> V=0.0

.....................................

-------- T = 1.0 ----------------- lower wall

The units are SI unit. Length=80, grid number in x=80, height = 1.0 and grid number is 20.

No-slip on upper and lower walls.

Pr=0.7, viscos=1.0e-4, Cp=1.0, Density=1. Converged steady upper half temp. field solution are listed here.

Solution method: SIMPLE, Co-allocated grid, 2-central difference scheme.

22 0.00E+00 1.00E+00 1.00E+00 1.00E+00 1.00E+00 1.00E+00 1.00E+00 1.00E+00 1.00E+00 1.00E+00

21 0.00E+00 4.44E-01 7.59E-01 8.32E-01 8.57E-01 8.72E-01 8.82E-01 8.90E-01 8.96E-01 9.01E-01

20 0.00E+00 1.18E-01 3.76E-01 5.06E-01 5.74E-01 6.17E-01 6.47E-01 6.70E-01 6.88E-01 7.04E-01

19 0.00E+00 1.29E-02 1.50E-01 2.57E-01 3.32E-01 3.88E-01 4.31E-01 4.64E-01 4.92E-01 5.15E-01

18 0.00E+00 -9.69E-03 4.51E-02 1.08E-01 1.65E-01 2.14E-01 2.56E-01 2.92E-01 3.22E-01 3.49E-01

17 0.00E+00 -7.84E-03 7.21E-03 3.58E-02 6.86E-02 1.02E-01 1.34E-01 1.64E-01 1.92E-01 2.17E-01

16 0.00E+00 -3.27E-03 -1.75E-03 8.08E-03 2.32E-02 4.17E-02 6.16E-02 8.25E-02 1.03E-01 1.23E-01

15 0.00E+00 -8.51E-04 -1.87E-03 3.16E-04 5.74E-03 1.42E-02 2.45E-02 3.67E-02 4.95E-02 6.30E-02

14 0.00E+00 -1.13E-04 -7.53E-04 -7.28E-04 6.02E-04 3.78E-03 8.12E-03 1.43E-02 2.12E-02 2.93E-02

13 0.00E+00 -2.24E-05 -1.28E-04 -3.65E-04 -3.36E-04 6.74E-04 2.05E-03 4.82E-03 8.10E-03 1.26E-02

12 0.00E+00 -5.39E-05 6.65E-05 -7.32E-05 -3.28E-04 7.77E-06 2.68E-04 1.62E-03 3.25E-03 5.94E-03

Hello, Mo Yun.

>The values at 4 corners are not used in FV method.

This does not matter. What matters is the behaviour of the derivatives of the solution when the corner is approched.

>Length=80, grid number in x=80, height = 1.0 and grid number is 20. Pr=0.7, viscos=1.0e-4, Cp=1.0, Density=1.

Hence, the grid Reynolds number is 1.0e+4? And central differencing. What about the stability?

Well, you see, it seems obvious that if I would continue explanations they were quite long. Reread textbooks on numerical methods. Pay attention to:

1. Central theorem of numerical analysis (though in your particular book it can be given another name or no name at all): Convergence follows from Approximation and Stability.

2. Note, that approximation depends on the properties of the solution: Namely, you consider quantities like (d^nU/dx^n)*h**2 as being second order. If, however, the derivative tends to infinity near a certain point, this is not so.

3. Find a place where grid Reynolds number is mentioned in the book and understand all about it.

Well, you see, there is no king's way in a science, as that old story tells. In a sense, this is why science is so attractive.

Good luck. And do not hesitate to ask questions if necessary.

Yours Sergei

(1). Before you formulate the problem, you need to know two things: one is the Reynolds number of your problem, and the other is whether is flow is going to be laminar or turbulent. (2). Since you were able to get some results, and don't want to get too complicated in using so-called advanced schemes, you can still run some additional cases and explore the results. (3). You can lower the Reynolds number by lowering either the velocity at the inlet of the viscosity. Try to run cases with Reynolds number=1, 10, 20,40,100,200 and so forth, and look at your results. (4). When you see negative temperature, try to increase the mesh density, say, from 20 to 40. And examine the difference between the two results. In 2-D problem, I think you should be able to increase the mesh density to a couple of hundreds without any difficulty. (5).Assuming that there is no error in your program, all of your results should be consistent and mathematically correct. ( consistent with the mesh used) (6). There is no problem with the negative temperature, the problem is with the mesh size and related truncation errors. at this point you may want to review the archived message posted some times ago under the title:CELL REYNOLDS NUMBER in this forum. (7). I am sure that your problem can be resolved with finer meshes, but ,if you don't want to invest in the finer mesh, the alternative is to use higher order up-wind methods ( or the like). (8). I don't think there is anything wrong with the results using the central differencing methods. Especially that you are not telling the equations to exclude the negative temperature values in the problem formulation.

Dear Sergei and Chien,

I deeply appreciate your suggestion. The Re should be defined as

Re=U_inlet * Height/ (density* viscosity).

Therefore, the Renolds number is 1000 in the case. I also checked the other Renolds numbers from 10 to 100. When Re is lower than 40, no negative results occurred. As we know, the laminar flow region under 2000 in the experiment knowedge, the numerical limit should be larger than this limit since the instable factor is smaller than the experiment.

As I mentioned, when I swithed the upwind scheme, that is the ** Temperature ** at the CV faces in the convection term Partial (u T)/ Partial (x) * Volume of cell are estimated from the upwind Cell center's value. The solution seems resonable but the numerical resolution is low theoretically. In my practice, CDS = (CDS + UDS - UDS)= UDS + (CDS-UDS), where CDS represents Central Difference Scheme and UDS does Upwind Difference Scheme. The steady and implicit scheme are used and the energy transport equation is discretized by UDS in the coefficients and (CDS-UDS) terms are dropped in the source term. Therefore, the solution is quite stable but remains the second accuracy.

As Sergei said, the combination of forward derivative and the backward derivative could be used for diffusion term discretization such as Partial(K Partial T/ Partial x) / Partial x etc but in the convetive terms, it denpends on the direction of velocity on the face of the CV. As for the diffusion terms, no problems occurs.

The effect of mesh resolution was also checked but same results obtained. In my experience, the solution should have no matter with the mesh resolution if we could put higher order accuracy bondary conditions in the diffusion equations. The analytical soltion is **just ** as the numerical one. But this rule does not suitable for convection.

In the energy equations, does the Cell Renolds number effect the temperature transport? In the explicit, transient method, the time step should be taken as small as the transport distance does not cross the cell because the linear approximation is valid only in the cell. But in this case, the flow field in the inlet that I got seems right, that is, the velocity concentrates to the central line of the duct and develops to the parabola shape in the lowere renolds number and thin boundary layers in the vicinity of the walls. In the discetized eq., the minus factor is in the (CDS - UDS) terms. In the upwind scheme, this term dispears. In the central diffence scheme, a heat sink is brought. Though the sourrouding temper is larger than zero but the nimus result might be given due to large heat sink. Of course we could multiply a blending factor f as f*(CDS-UDS) where f=0 to 1. But this method does not stand on the theoretical ground, only the experience. Maybe the probelm is in how we should formulate the blending factor f theretically, so that no nonphysical results occurs and with higher order occuracy.

Sinerely yours,

Yun Mo

(1) Based on Re=40, one can compute the cell Reynolds number by dividing the total number of cells across the channel,that is Re,c=40/20=2. (assuming that the cells are equally spaced) This Re,c=2 is a magic number for the central difference scheme when used for the convection term. So, I would say that the negative values are related to the so-called cell Reynolds number limit when using central difference for the convection term. (2). Regardless of the methods used, (you are free to create any new scheme you like)as long as the solution is grid independent , the solution is all right. Even with the up-wind method, you should be able to get improved accuracy by reducing the mesh size.

Dear Chien,

I would very appreciate it that you have give me the idea. Yes, Cell Reynolds Number really decides the transported quantity when the central diffence scheme is used. According some textbooks, it also called Peclet number. Therefore, the cell number is defined as

Re, cell = U* dx/ nu

where dx is transport length in the cell in U direction,

nu is viscosity divided by the density,

but not the one of Duct Reynolds number divided by the cell number in the cross section. I reduced the length from 40 to 1 in the flow direction and remained the number of division in flow direction. The temperature field **really** changed from minus values to plus values around the cell reynold number Re=u * DX/ NU =2. Therefore, all of schemes such as upwind scheme, Hybrid scheme, exponential scheme and QUICK method are the praticed method to deal with this kind of probelm. Since then, when in the multiple grid method, the best way to start from the upwind wind (or hybrid method) and the central difference scheme can only be switched when the Cell Reynolds number becomes less than 2.

Thank you all very much. This problems may be concluded here.