Hello,

I have the following equation:

C1 dU/dt + v d/dx( C2*U) = d/dx(C3*dU/dx)

U= continuous function v= constant (let's say ~0.008) C1, C2, & C3 are variable coefficients with discontinuity.

I am solving this equation in a 1-D domain. I am expecting 3 different regions inside this domain (phase-change flow problem). In the 1st and 3rd regions, C1, C2 & C3 are all constants. In the 2nd region, only C2 & C3 change with position and undergo abrubt & discontinuities changes at the transition locations. I tried many methods (e.g. central differencing, upwind differencing, & MacCormack methods) with no luck. Could you suggest a method for solving such a problem? Do you expect your method will work for a 2-D problem?

More info: 1st region: 6.5e9<U<6.0e9

C1~1.0, C2~0.001, C3~1e-7 2nd region: 6.0e9<U<6.0e8

C1~0.4

0.003<C2<4.0 4.0e-5<C3<0.0 3rd region: 6.0e8<U<1.0e6

C1~2500, C2~1.6, C3~2.0e-4

Many Thanks!

This appears to be a very stiff problem. Assuming you set the system of equations properly, when you say you tried different methods with no luck, what do you mean? Are you getting stability problems or else? If the former, have you looked at the stability conditions (analytically) or experimented with smaller and smaller time and spatial scales and check for potential trends?

Adrin Gharakhani

Dear Adrin Gharakhani, Many thanks for your reply. What I am trying to solve is only a single equation with its proper boundary conditions (Dirichlet & mixed). I have found that for a small time step (1e-8) {this is experimental check and it is very hard to study stability conditions with such variable coefficients and multiple discontinuities}, the problem works fine but it takes a lot of time to converge, obviously. If the time step is greater than "1e-8" I will get oscilatory (and overshoot)solution with no physical meaning which is due to instability condition(s) "as I think". I am not sure if the mixed B.C. [a*dU/dx + b*U = c] has to be handled carefully especially with variable coefficients. However, this small time step can be increased to (0.0001) in the other two regions but I have to maintain the small time step in one region.

Regards

PS I have no idea about "check for potential trends"

> PS I have no idea about "check for potential trends"

Trends such as conversion vs. oscillation using different timesteps, etc. You have already answered your own question. There is a timestep size above which you get oscillatory solution. If you were to combine that timestep with the local grid spacing and characteristic velocity (or diffusivity, etc.) you can see some type of a stability condition. The fact that you get a solution with small timesteps is positive and suggests that the methods that you are using have this inherent limitation for stability. So, really, you have little to be concerned with. If you want larger timesteps you need to look into the literature and find ways for solving highly stiff problems. Mixed boundary conditions are definitely more difficult to converge. Just make sure that you have the correct match across different regions. One suggestion, you can use adaptive timestepping based on some local stability condition, as well as "subcycling" - for regions where you can use large timestep, "freeze" the time advancement temporarily and march the region requiring more timesteps in time. For example for region 1 the timestep may be dt1 and for region 2 the timestep may be dt2. Then for each dt1 you have to use dt1/dt2 steps to move the stiff part of the solution in time.

Adrin Gharakhani

If you you want to capture your discontinuity with fair accuracy and without the oscillations you will need to move to a high resolution scheme.

I would suggest reading: randall, LeVeque Numerical Methods in Conservation Laws...

(Or search on the Web for his class notes...)

(1). Since the equation is not valid at the interface where the coefficient is discontinuous, it is a good idea to solve three problems instead of just one. (2). If you expand the equation, then you can see easily the singularities associated with the derivatives of the coefficients C2, and C3. (namely dC2/dx, and dC3/dx)

Dear Prof. John C. Chien,

The equation is valid throughout the different regions "including" interfaces and discontinous points.

Regards

(1). Well, if that is the case, then you should have no problem at all.

Dear Prof. John C. Chien,

Obviously, I mis-represent what it is my real problem. Although the equation is valid throughout the different regions, the coefficients are highly non-linear and nead to be handled carefully (which I don't know how). However, there is an elegant suggestion by Mr. Gharakhani that I am trying now and I hope it works for my case. What I really like to know if the Upwind scheme or MacCormack scheme can handle this or not. Or, if anyone could suggest another (finite difference FD) scheme.

Regards