[hunt_mat]

Originally Posted by FMDenaro

Your system can be written like


u^n+1 = T.u^n

Where u is the state vector.

If you do the stability analysis on the eigenvalues of T you will be able to understand that the CFL is not the only parameter to check.

I'm seeing that. I'll do the stability analysis to see if I can't get an answer.

[LuckyTran]

Sure, it is a system of equations with a trivial method to solve it. I just wanted to make sure whether you have a system like Ax=B (this to me is a linear system of equations) or a system like X=C (this to me is an assignment operation).

I don't know elastic beams well enough to know whether your Peclet number is greater or less than 1. But for high Peclet numbers I don't recommend forward differencing on advection terms (this is downwinding). I would recommend central differencing or backward differencing (backwinding) unless you have Peclet number around 0.

Central differencing on the diffusion terms is reasonable. As hinted earlier, for Peclet numbers not greater than 1, diffusion terms are unstable with Fourier number.

Have you considered debugging on a trivial case? I.e. homogeneous BCs, zero U, and showing that you recover the expected solution and that there is no silly typo? This test case would also be helpful since you can arbitraily coarsen your grid to make the problem more stable to find such bugs.

[FMDenaro]

Originally Posted by hunt_mat

I'm seeing that. I'll do the stability analysis to see if I can't get an answer.

Could you show the analytical solution you have?

[hunt_mat]

I am now examining stability conditions using the Lax equivalence theorem and I am getting some difficult conditions that don't seem to be consistent which is very worrying.

I've included the beginning of my analysis.

[hunt_mat]

Originally Posted by LuckyTran

Sure, it is a system of equations with a trivial method to solve it. I just wanted to make sure whether you have a system like Ax=B (this to me is a linear system of equations) or a system like X=C (this to me is an assignment operation).

I don't know elastic beams well enough to know whether your Peclet number is greater or less than 1. But for high Peclet numbers I don't recommend forward differencing on advection terms (this is downwinding). I would recommend central differencing or backward differencing (backwinding) unless you have Peclet number around 0.

Central differencing on the diffusion terms is reasonable. As hinted earlier, for Peclet numbers not greater than 1, diffusion terms are unstable with Fourier number.

Have you considered debugging on a trivial case? I.e. homogeneous BCs, zero U, and showing that you recover the expected solution and that there is no silly typo? This test case would also be helpful since you can arbitraily coarsen your grid to make the problem more stable to find such bugs.

An operator A is called linear if and only if A(au+bv)=aAu+bAv, the operators in my original equations satisfy this condition.

I'm currently trying to derive such conditions through a stability analysis using the Lax equivalence theorem. I do not have an analytical solution, it would have to be through a Fourier series. You can't do a separation of variables in this case as the system is coupled. I think that even looking for an analytical solution may be difficult. I could look for the case for a constant temperature, this would still be physically realistic and the analytical solution might be easier to obtain.

[hunt_mat]

Originally Posted by FMDenaro

Could you show the analytical solution you have?

I don't have one.

[LuckyTran]

Look man, we all know this is a linearized form of a more complex set of governing equations. I would hope that your linearized system is indeed linear. From the start I've made it clear that I was asking whether you have an implicit system of algebraic equations that need to be solved using gaussian elimination or if you have an explicit time-marching scheme. Anyway, thank you for confirming that you have both.


a5/a4 is your Peclet number

[FMDenaro]

Originally Posted by hunt_mat

An operator A is called linear if and only if A(au+bv)=aAu+bAv, the operators in my original equations satisfy this condition.

I'm currently trying to derive such conditions through a stability analysis using the Lax equivalence theorem. I do not have an analytical solution, it would have to be through a Fourier series. You can't do a separation of variables in this case as the system is coupled. I think that even looking for an analytical solution may be difficult. I could look for the case for a constant temperature, this would still be physically realistic and the analytical solution might be easier to obtain.

But have you dobe first the stability analsyis for the case of a single equation of convection-diffusion-production ?

[hunt_mat]

Originally Posted by LuckyTran

Look man, we all know this is a linearized form of a more complex set of governing equations. I would hope that your linearized system is indeed linear. From the start I've made it clear that I was asking whether you have an implicit system of algebraic equations that need to be solved using gaussian elimination or if you have an explicit time-marching scheme. Anyway, thank you for confirming that you have both.


a5/a4 is your Peclet number

It wasn't linearised from any set of equations, I simply wrote it as a linear system that has all the problems of the nonlinear system that I was working on. I have looked at this from several perspectives: 1) Using the backward Euler method, this should be unconditionally stable, it wasn't. So then I did the simple solution, forward Euler method, I have such stencils in the pdf I attached in a previous post. Nothing seems to be working. I have done a stability analysis using the norm of the matrix being less than 1, and I have gotten some problems with this as I cannot get a consistent set of conditions.

I'm begnning to lean on my original was inconsistent.

[hunt_mat]

Originally Posted by FMDenaro

But have you dobe first the stability analsyis for the case of a single equation of convection-diffusion-production ?

I did do a stability analysis, I found that if c=Udt/dx, and a=Ddt/dx^2, then I have the following conditions:
a>c and a<1/2 gave a stable solution along with c>a and a<1.

I found this out by examining the norm of the matrix and setting it to be less than 1 for stability.

[FMDenaro]

Originally Posted by hunt_mat

I did do a stability analysis, I found that if c=Udt/dx, and a=Ddt/dx^2, then I have the following conditions:
a>c and a<1/2 gave a stable solution along with c>a and a<1.

I found this out by examining the norm of the matrix and setting it to be less than 1 for stability.

No, first you have to consider also the production term, then is better you do then von Neumann analysis and evaluate the full stability regions (it is an isosurface T |G|=1 in a 3D space.

[hunt_mat]

Originally Posted by FMDenaro

No, first you have to consider also the production term, then is better you do then von Neumann analysis and evaluate the full stability regions (it is an isosurface T |G|=1 in a 3D space.

Why do you think using the Lax equivalency theorem is wrong in this situation?

[FMDenaro]

Originally Posted by hunt_mat

Why do you think using the Lax equivalency theorem is wrong in this situation?

Lax equivalent theorem is valid. The issue is that you have to prove the stability.
Consider that single equation, you have

df/dt+u*df/fx=d*d^2 f/dx2 + k*f

The stability region is obtained as the isosurface:

|G|(cfl, Re_h, gamma)=1

G being the amplification factor of the scheme.

That is you have a functional relation between the various parameters.

[hunt_mat]

Originally Posted by FMDenaro

Lax equivalent theorem is valid. The issue is that you have to prove the stability.
Consider that single equation, you have

df/dt+u*df/fx=d*d^2 f/dx2 + k*f

The stability region is obtained as the isosurface:

|G|(cfl, Re_h, gamma)=1

G being the amplification factor of the scheme.

That is you have a functional relation between the various parameters.

von Neumann stability analysis isn't always the easiest thing to do.

[FMDenaro]

Originally Posted by hunt_mat

von Neumann stability analysis isn't always the easiest thing to do.

Are you kidding?? It is the MUST in any numerical study! Or is a fundamental topic in a CFD course.

For a single 1d equation is very very easy to do using matlab or maple software, just few lines of code.
See https://www.researchgate.net/publica...tion_Equations

[hunt_mat]

Originally Posted by FMDenaro

Are you kidding?? It is the MUST in any numerical study! Or is a fundamental topic in a CFD course.

For a single 1d equation is very very easy to do using matlab or maple software, just few lines of code.
See https://www.researchgate.net/publica...tion_Equations

Code? All the stability analysis I've seen has been done algebraically.

[FMDenaro]

Originally Posted by hunt_mat

Code? All the stability analysis I've seen has been done algebraically.

In order to generate the practical 3d map for the 1d equation, you have to investigate the most critical frequency for each value of CFL, Re_h and gamma.
Indeed, the complex amplification factor depends on 4 independent variables. Are you able to generate that algebrically ? Let us see …

[hunt_mat]

Originally Posted by FMDenaro

In order to generate the practical 3d map for the 1d equation, you have to investigate the most critical frequency for each value of CFL, Re_h and gamma.
Indeed, the complex amplification factor depends on 4 independent variables. Are you able to generate that algebrically ? Let us see …

In the 1D cases I've seen, the von Neuman stability analysis reveals an inequality on the parameters involved in the stencil. I don't know what you're talking about, the type of stability analysis I'm using is setting v_{n,m}=e^{a*n*dt}e^{i*b*m*dx} and looking for the condition for |e^{a*dt}|^{2}<1.

[FMDenaro]

Originally Posted by hunt_mat

In the 1D cases I've seen, the von Neuman stability analysis reveals an inequality on the parameters involved in the stencil. I don't know what you're talking about, the type of stability analysis I'm using is setting v_{n,m}=e^{a*n*dt}e^{i*b*m*dx} and looking for the condition for |e^{a*dt}|^{2}<1.

No, not this way. Let me consider a simple FTUS scheme




f(i,n+1)=f(i,n)-cfl*[f(i,n)-f(i-1,n)]+(cfl/Re_h)*[f(i+1,n)-2*f(i,n)+f(i-1,n)]+gamma*f(i,n)


cfl=u*dt/h, Re_h=u*h/ni, gamma=dt*k



That is



f(i,n+1)=[1-cfl-2*(cfl/Re_h)+gamma]*f(i,n) +[cfl+(cfl/Re_h)]f(i-1,n)+(cfl/Re_h)*f(i+1,n)


Hence


G=[1-cfl-2*(cfl/Re_h)+gamma]+[cfl+(cfl/Re_h)]*exp(-i*k*h)+(cfl/Re_h)*exp(-i*k*h)




From this expression you evaluate the expression for the modulus |G|(cfl,Re_h,gamma,k*h)


so that you can simply code the research for all -pi<=k*h<=+pi for the most critical value of |G| at values clf,Re_h,gamma.


That can be easily done by saving the max value of |G| in a 3d matrix and then extracting the isosurface at |G|max=1 value. Slices can help to see the 2D maps.


Very easy to code.

[hunt_mat]

Originally Posted by FMDenaro

No, not this way. Let me consider a simple FTUS scheme




f(i,n+1)=f(i,n)-cfl*[f(i,n)-f(i-1,n)]+(cfl/Re_h)*[f(i+1,n)-2*f(i,n)+f(i-1,n)]+gamma*f(i,n)


cfl=u*dt/h, Re_h=u*h/ni, gamma=dt*k



That is



f(i,n+1)=[1-cfl-2*(cfl/Re_h)+gamma]*f(i,n) +[cfl+(cfl/Re_h)]f(i-1,n)+(cfl/Re_h)*f(i+1,n)


Hence


G=[1-cfl-2*(cfl/Re_h)+gamma]+[cfl+(cfl/Re_h)]*exp(-i*k*h)+(cfl/Re_h)*exp(-i*k*h)




From this expression you evaluate the expression for the modulus |G|(cfl,Re_h,gamma,k*h)


so that you can simply code the research for all -pi<=k*h<=+pi for the most critical value of |G| at values clf,Re_h,gamma.


That can be easily done by saving the max value of |G| in a 3d matrix and then extracting the isosurface at |G|max=1 value. Slices can help to see the 2D maps.


Very easy to code.

Then I think we're talking about different things. Could you give me a link regarding the stability you're talking about?