it must be clearly defined in numerical analysis even it is right in purely mathematical theory. How do you expect the numerical solution for the first time step if you put singularity condition in it? Teach the machine a way.

You need to understand the domain of validity. In your problem, your can either say that (i) you have delta function of BC: T(0,t) = delta(t) where detla(t) means at t=0 it is zero and at t>0 it jumps to 1, or (ii) you have an IC that is not valid at the boundary (x=0). In either way, there is no problem in the analytical solutions.

However, when you work on numerical solutions, you need a BC and an IC to be valid at the SAME TIME, say at t=0. I don't know exactly how you solved the problem, but if you apply constant BC at x=0 as T(0,t)=1 from t=0 and up, you are not solving the same problem as your previously defined. You are solving something like

IC. T(x,0)=1, x=0; T(x,0)=0, x>0.

which is inconsistent with your IC for the analytical solution.

Or, if you use time-variable B.C. you are solving something like

B.C. T(0,t)=0, t=0; T(0,t)=1, t>0

which is different form your B.C. in the analytical solution.

It is actually very easy to numerically solve such a problem if you correctly setup the BC or IC. Depends on your understanding of the difference of t=0 and t>0, you may have difficulties to get the corresponding analytical solution.

Many real world problems involve singularity for fluxes and tractions, so there is nothing wrong in numerically solving these singularity-type problems. Think of shock waves, as concerns CFD. Or, even easier, consider a stress singularity right at the step corner for the backward-facing step flow. Many CFD people solve this "benchmark problem" without realizing that, oh, wait a minute, do I have a singularity here?

We don't want to mess up with any other problem. Fine, when no singularity is present -- conventional numerical approaches do okay. But, staying with the problem initiated this thread, it involves something that number-crunching people hate.

So please, simply state -- can you solve or not? I'm wasting my time otherwise.

Sorry, forgot to mention that it's not a delta function involved here, rather a Heaviside step function, so that somebody would not be misled.

(1). I think the problem is actually a model. An idealized model of the limiting process. (2). In physics, it is not possible to have the state change in zero time, and it is not possible to have a state jump at the same point in space. (3). Since a point in space has no dimension(or size), the point x=0, and the point x=0+ is basically the same point, unless a finite value is given. For example, x=0 + 0.0000000001 . (4). This is also true for the time variable. time=0 , and time=0+ is the same time. Mathematically, it is also not possible. But under the limiting processes, it can be defined mathematically. That is to represent the limiting case. (5). As a consequence, the solution also must be expressed in terms of the limiting process.

Don't mess up the definition of singularity for different situations. You are giving unconsistant B.C. and I.C. to start a numerical analysis problem which require a constant B.C. and I.C. at a cross point. It is nothing to do with back step and shock wave. You are also interested in the first time step. You started with unconsistance and look for a consistant result out of it.

thanx Jonh. You are right.

Having been involved with "traditional" (not advanced Mikhail ) and advanced CFD for a while, I can clearly understand why people responding to this thread would simply dismiss the problem as incorrect/incosistent without giving it some serious thought first.

Most traditional CFD techniques shove these singularity issues (such as those due to sharp corners) under the rug! This is usually justified by noting that in real life there is no singularity (sharp edges). However, if there is no sharp edge, then there is some curvature that the control volumes that are used in CFD (with their own sharp edges) do not emulate properly. Thus, you can never honestly say what exact problem it is you are solving with this type of an approach. As a matter of fact, if we were to conduct a "serious" test for grid dependence, soon those numerical techniques that do not handle singularities properly will break down as the grids get smaller and smaller (and we get closer and closer to the sharp edges). Don't expect the solution to get better as the grid density increases significantly near the sharp edges! What one may observe is a leveling off of the solution because it hits the precision limit of the computer, but that does not imply convergence!

Now, back to the main question. I agree that there is no inconsistency in the posed problem - there is just a discontinuity in one of the BC's. I think there would be less confusion had Mikhail posed the problem using the Heaviside step function argument.

T(0,t) = H(t) t>=0 T(1,t) = 0 t>=0 H(0-)=0 H(0+)=1

This is definitely a wellposed problem.

It's like having a rod with zero temperature, and all of a sudden, in such a short time that up to the machine precision dt=0, applying temperature to one end.

The argument that in real life there is some time lag is not a good one. Numerical techniques don't know whether a well-posed problem is physical or not. They are expected to faithfully get the discrete (digital) solution to a given analytical (analog) differential equation. Period.

The original question was whether people on this list are aware of such a method that can get the discrete solution to the stated problem with discontinuous BC.

Mikhail, I understand you cannot divulge the information right now, but I would like to learn how you handle the jump (or singularity) when you can.

Thanks for the good question

Adrin Gharakhani

Dr. Gharakhani:

Thank you for your input, I really appreciate it as well as many people here will.

As for advanced methods, I should have used an apostrophe in the word "advanced" -- that was nothing but an irony.

Mikhail

If you impose the BC and IC that way (I call it the revised problem), you will not get the same analytical solution! I totally agree with you that your BC/IC are well posted, and the problem is reduced to a discontinity issue. I remember the question was how to get a first time step solution that has less than 10% error when compared with the analytical solution with the original BC/IC. If you change the BC/IC, you are working on a totally different problem and you will get different analytical solution. That is the whole point of inconstency.

In the original set of BC/IC, there is no way to solve it numerically, because at t=0, the temperature field was not defined at boundary x=0. In order to solve the problem, you MUST specify a value of T(0,0), but when you do that, you changed the problem specification. This can be easily verified by checking the analytical solution using x=0 and t=0, do you get the same answer as you specified? Of course not. So the error does not come from your first time step, it comes even before you started it. In such a problem, there is no consistent numerical solutions for the corresponding analytical solution.

In the revised problem with Heaviside (this is my mistake, thanks.) IC or BC, the problem is a discontinuity issue. Mathematically t=0 and t=0+ is unable to numerically model If we understand that numerical modelling is an approximation process, we can alwasy model t=0+ as a very small t0 and then devide this periodic by even small dt. It is something like we model shock waves in space domain. While we know that shock wave thickness is the order of 1/SQRT(RE), our best approach is at leat one grid point thickness.

Back to the main story. It is a very good idea to compare numerical solutions with analytical solution to check the basics of numerical scheme. But we must make sure we are comparing an apple to an apple.

You have led people misunderstanding the problem. You just gave a description of the problem but have not told us how did you solve the problem.

Let me make it clear.

This is a transient diffusion problem, dT/dt=kd^2T/dx^2, with boundary condition: x=0, T(0,t)=1; x=1, T(1,t)=0; and initial condition: T(x,0)=0.

I think you have not note the fact that the initial condition is singular, and its first and second spatial derivatives do not exist at all. If you had used a explicit scheme to descrite the equation, you could not have got a correct solution certainly, because you can not calculate the spatial derivatives at the boundary x=0 where dT/dx=infinite on t=0. Fortunately, the solution of t>0 can be analytically obtained and then we know that the spatial derivatives must exist when t>0. So for this impulsively started problem, a fully implicit scheme should be used!!! I guess that you must have calculated the spatial derivatives on the instant t=0. If so, this is really wrong. Your computation is absolutely unstable. You'd better read some textbook about the stability of explicit scheme.

If you tried to calculate this problem by using the following implicit formulation, you would find that the world is beautiful.

Tnew - k(dt/dx**2)(Tnew(i-1)-2Tnew(i)+Tnew(i+1)) = Told

As I can understand, the solution must be a time-developping boundary layer. The temperature is developped as the time marches towards and becomes steady when t -> infinite. You should use a small time step and small spatial size to obtain a accurate solution.

One should simultaneously consider the consistency of the descrite scheme, boundary condition, initial condition and so on when he solve a problem. In fact, many commercial softwares use explicit scheme to solve such application problems. As we know in this discussion, this is really dangerous. Many numerical shemes are not suitable to problems with singularities of boundary condition and initial condition.

I believe that at least for this simplest problem, numerical computation (or CFD here) can get a correct solution if a correct scheme is used.

Hope this help you for understanding the problem.

Z. Lei

John, I can give you a dozen similar problems like this.

For example, the flow near the sharp leading-edge of a flat, if you take t -> x, and x -> y here, you will find this equation becomes a thin boundary layer equation.

Another example is stagnation flow near a wall with a constant temperature on t>=0 and zero in the inner field. The problem may be two or three dimensions.

And then many impulsively started problems, including shock wave problems.

Read another post of mine, you will know why he can not get a correct solution and why he is wrong.

Z. Lei

(1). I think, the basic problem with this transient diffusion heat conduction problem is that at time=0, the solution is not unique. (2). Mathematically, the condition can be defined from (time=0-) , (x=0-) and (time=0+), (x=0+), only. (3). In the case of a shock, the upstream condition must be supersonic (Mach>1), and the downstream solution must be subsonic (Mach<1). The solution at the shock therefore can only be approached from both side, but not at the location of the shock. (One could solve for the continuous shock structure problem though). (4). With this limitation, there is no way one can include the time=0 and x=0 as part of the solution. The heat flux will be either zero or infinite. (5). The reason why the solution exist is because this is a diffusion problem . The solution can therefore be computed from the time=0++ side. If it is a convection problem, one must use something like method of characteristics to follow the discontinuity, rather than across the discontinuity. (6). The other way to do is to replace the discontinuity by an infinite series of continuous functions. Perhaps this would be a more meaningful boundary condition. (7). Using the implicit method is also a good approach in the diffusion problem. The diffused solution would always stabilize the calculation. (8). As for the use of the explicit method, the most commom troublesome problem is the initial start of a flow field solution. Because of the poor initial guess, and the use of extremely small mesh near the wall with a low Reynolds number model, the initial gradient can be extremely large. In most cases, the use of explicit method will give diverged solution within one iteration or time step. This is a real problem. (9). My feeling is that, if there is a solution to this problem, numerical or analytical, the solution must come from the time=0++ side because one can't march through the time=0 singularity.(I though I said somewhere before, it is a black point. In a sense that one can't march through a black hole. And modern calculus is basically mathematics of limiting processes. At time=0, the H-function is not unique.)

John, I agree with your opinions. As metioned in my last post, in fact the intitial condition does not satisfied the control equation. One can not discrete the equation on the instant t=0. The discrete should be started from t=0+. Whatever the discrete scheme is used, BEM, FDM, FVD and FEM, the computation should carry out in this way,

(Tnew - Told)/dt = K (d^2T/dx^2)new

to avoid the singularity of the initial condition.

This is really a very simple question, and can be properly calculated without any problem. It has been discussed as an example in many textbooks. And, similar problems have been also solved in many practical applications.

Z. Lei

If your erudition doesn't go beyond finite differences, that's your problem -- definitely, not mine. Sure thing, all this garbage-makers like FDM, FVM have no way to take care of singularity of the problem, so you simply ignore it. That's fine with me. But ignoring the problem is not the best idea at all...

Period.

(1). I would say that solving a problem with singularity is largely outside the range of traditional CFD methods, which must deal with well-defined conditions and continuous solutions. (2). The classical series expansion method of solution or the singular perturbation type of analysis, does not enter the CFD formulation. (3). The problem is real with a Laser tool. Since the time is so short , the temperature gradient is so steep, and the heat flux is so large, before the material can dissipate the heat by conduction, it simply evaporate.

i find this mikhail character to be extremely rude. he poses a physically impossible problem then he refuses to listen to mr Zhong Lei when he proposes a technique to solve his problem. hopefully this man is not an engineer after all we engineers typically do not concern ourselves with non-physical problems unless they are a simplification. i'll trust mr. Lei who say's he can solve the problem. but whether he can or not is irrelevant because real substances cannot tolerate infinite temperature gradients or rates of change of temperature and neither can most numerical techniques. if he's hot to solve this problem use the analytical technique. it's the right tool for the job. as i have said before i'm not a fan of solving laminar, incompressible pipe flow with 4 stage runge kutta. he should use the appropriate tool for this silly problem, wgich is the pencil and paper not the computer. and if people with more skill than himself like Mr. Lei answer his question he should listen.

Clifford, I would tend to agree with you.

One question that nobody sems to have asked yet : what is the point in getting an analytically accurate solution of the problem at the first time step ??? Other than that of the brief intellectual satisfaction one might get.

Yah, like a fussy baby. Refused logic reality and attacked helpers.