# Why cant I get a stable Cranck Nicolson discretization?

2 visualizzazioni (ultimi 30 giorni)
Jort Puiman il 26 Gen 2023
Commentato: Jort Puiman il 22 Mar 2023
Hi all,
For a project we try to design a reactor with a diffusion-convection-reaction mechanism before it reaches a steady state. Now, we are using the Cranck Nicoloson discretization for this reaction. However this does not seem to work properly. We get an unstable system, because we have a low diffusion coefficient. We have been playing around, but we cannot find how we can make this system more stable.
Attached is the code we are using right now.
Do you maybe have suggestions on how to make the discretization more stable? Many thanks!
edit: the initial mass balance we are trying to solve is dx/dt = D*d^2x/da^2+v_av*void*dx/da-k*void*x
We use a Dirichet boundary condition so that x=0 for t = 0 and a = 0.
##### 4 CommentiMostra 2 commenti meno recentiNascondi 2 commenti meno recenti
Bjorn Gustavsson il 21 Mar 2023
@Jort Puiman: were our suggestions sufficient to solve your issue?
Jort Puiman il 22 Mar 2023
Unfortunately, our model was too sensitive, and our lecturer made us switch to the method of lines discretization.
The timesteps were indeed problematic, but smaller timesteps did not fully fix the problem, and the model remained unstable after a longer time.

Accedi per commentare.

### Risposte (2)

Bjorn Gustavsson il 26 Gen 2023
(caveat: I didn't run your code or analyze it, but speak from my experience with CN) Most likely you take too long steps in time. Even for CN-methods one has to obey the CFL-condition when selecting the time-steps. You might have to read up on the Courant-Friedrichs-Levy condition: wiki-on-CFL-condition.
That was my experience, didn't think I needed to bother because CN is stable, but first solutions were not good, when the CFL-condition on dt was adhered to all was good.
HTH
##### 0 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

Torsten il 26 Gen 2023
Spostato: Torsten il 26 Gen 2023
We use a Dirichet boundary condition so that x=0 for t = 0 and a = 0.
And which conditions at a = a_end ?
Try "pdepe".
If the convection term v_av*void*dx/da is too strong compared to the diffusion term, discretize dx/da using an upwind method, not central differences.
##### 0 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Categorie

Scopri di più su Geometry and Mesh in Help Center e File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by