
#1
Feb2112, 04:49 PM

P: 9

I have a perfectly working 2D finitedifference thermal solver, using an alternatingdirection implicit scheme (painfully written in Fortran 90 at my adviser’s insistence), that I have recently extended to 3D. It is essentially a three step approximation to the CrankNicolson equation, which is shown to be unconditionally stable in homogeneous media with constant boundary conditions in the old text book in which I found it.
Who the hell writes a thermal solver to simulate a single material with constant temperatures at the boundaries? Anyway, I have adapted the method to an inhomogeneous material structure, and, regardless of boundary conditions, I have exponential error growth (in magnitude...the actual error oscillates), which starts at the material interfaces. For almost prohibitively low timesteps, though still much larger than what would be required for an explicit method, the thermal solver actually appears stable and produces the correct results as far as I can tell, so I don't think there is any simple error in my implementation. It appears to me that the material inhomogeneities have broken the unconditional stability of the 3D alternating direction implicit scheme, and I am having a very hard time finding relevant resources. Posting this here is probably a bit of a fishing expedition, but if this sounds familiar to anyone else, please elaborate on your experiences. Thanks in advance. 



#2
Feb2112, 06:59 PM

Engineering
Sci Advisor
HW Helper
Thanks
P: 6,348

The first thing to check is if it works in 3D with a homogeneous material.
The "point" of ADI is that you are taking steps in different directiions which are unstable, except (at least in 2D) the instabilities in the two directions cancel out. I don't know how that would extend either to 3D or inhomogeneous materials, but (leaving aside bugs in the code) I guess that's one place to start researching. FWIW explicit forwarddifference finite difference schemes (with low time steps for accuracy) were used for 3D problems back in the 1970s, but things have moved on a bit since then. For an explicit scheme, you can do the calculations pointbypoint, so you don't need anything like "ADI" to avoid solving sets of linear equations. You can generalize Crank Nicholson into a family of methods with a parameter (say ##\theta##) where ##\theta = 0## is forward difference, ##\theta = 1/2## is central difference, and ##\theta = 1## is backward difference. Sometimes the best compromise of accuracy against stability is about ##\theta = 0.6## or ##\theta = 2/3##. (It's possible to invent an armwaving argument that says 2/3 is a "good" value to choose, but I'm not personally very convinced by it.) 



#3
Feb2112, 08:17 PM

P: 9

It works just fine with homogeneous media and simple boundary conditions, but that also would bypass a lot of potential errors in the code, e.g. typos in interpolating the staggered material mesh onto the temperature mesh.
For the system I need to simulate, which contains both glass and silver, and has 20 micron cylindrical substructures, I would need time steps on the order of femtoseconds for an explicit algorithm. It is starting to look like I can get away with 10's of nanoseconds with the ADI, but it certainly isn't unconditionally stable. I have a fully functioning inhomogeneous θ=1 1D straight CN, and a 2D ADI code that are stable in every scenario in which I have tested them, even coupled to radiative boundary conditions with source inhomogeneities (i.e, well outside the formal proof of unconditional stability), so I am not a complete novice on CrankNicolson implementations. There is just not much literature on real three dimensional ADI implementations. While the source for my algorithm (I can post it and a copy of the relevant paragraph tomorrow) claims it is unconditionally stable for the example given for homogeneous media, I am not familiar enough with numerical analysis to determine easily if the structural inhomogeneities should break the unconditional stability, as they appear to. Then again, maybe nobody is skilled enough to do it "easily." If in principal ADI's are stable due to canceling errors in each explicit part of the update, perhaps it is not general inhomogeneity that breaks the stability, but anisotropy. Some fully symmetric 3D test cases do seem to require a less refined timestep than the extrusion style test cases. 



#4
Feb2112, 09:12 PM

Engineering
Sci Advisor
HW Helper
Thanks
P: 6,348

CrankNicolson instability in 3D
Googling around, the DouglasGunn method seems to be the favorite 3D ADI method.
Re your test problems and interpolation etc, some problems where the steady state is a uniform heat flux in different directions, (including heat flux not parallel to the finite difference grid) might flush out some errors. Most of my heat transfer work has been with FE codes rather than FD, so unfortunately I've forgottem most of what I used to know about FD algorithms 



#5
Feb2512, 02:08 PM

P: 8

The 'fix' is the PeacemanRachford method. BTW do you have mixed derivatives like u_xy etc.? A better approach imo is the Method of Fractional Steps (Marchuk/Janenko). And the Alternating Direction Explicit (ADE) is used as well. Google will help. 



#6
May1812, 07:02 AM

P: 1

1. I suppose that the proof of the unconditional stability of the scheme that you mention is based on the von Neuman analysis... It is important to understand that this analysis does NOT take into account any boundary conditions! It gives the stability criteria only for the FD scheme itself. 2. I know that the generalizations of a 2D ADI scheme to 3D case, known as the PeacemanRachford scheme, is NOT unconditionally stable. One must consider the schemes known as DouglasRachford or BrianDouglas  these schemes are unconditionaly stable in 3D ADI. 3. Concerning the boundary conditions in ADI methods: there is a problem with the socalled intermediate boundary conditions  in a correct implementation one must specify boundary conditions at the intermediate time steps. For details see the book by Strikwerda "Finite difference schemes and partial differential equations"  there is a discussion of this problem for the 2D ADI scheme. The same is valid for 3D. The illposed boundary conditions can, in principle, produce instabilities. But if you have a nonlinear heat equation (e.g., thermal conductivity depends on temperature) the situation with imtermediate boundary conditions becomes much more complicated... 4. An anticipated question: is there a method of stability analysis that could take into account the boundary conditions? I do not know. But I suppose that the method of energy inequalities could solve this problem. Unfortunately I am not familiar with this method and I would be very grateful if someone could explain it to me. Probably one could apply a Matrix stability analysis. Please, let me know if you have advanced in solving your problem with inhomogeneties. Actually, I work on a similar problem  I have to pose 3rd order boundary conditions in some internal points of the modelling region, i.e. I also have inhomogeneties. good luck! 


Register to reply 
Related Discussions  
Is there more than one CrankNicolson scheme?  Differential Equations  1  
Crank Nicolson method  Programming & Computer Science  0  
Crank Nicolson Stability  Differential Equations  0  
CrankNicolson method (Matlab)  Math & Science Software  0  
Parallelizing CrankNicolson method  Differential Equations  3 