Finite volume methods

FMDenaro · August 10, 2023, 09:33

Quote:

Originally Posted by hunt_mat

The position starts with u=0, the Neumann condition is a stress free condition basically, it's the temperature that drives the process. I'm not too sure for nu, as the boundary condition comes from the velocity.

Ho do you express nu as a function of u? That seems more a sort of Robin condition, not exactly Neumann.

FMDenaro · August 10, 2023, 09:45

In order to get an idea of the Reynolds number Vr*L/D, one needs to estimate Vr from the BCs:

du/dx=-k*nu -> du*/dx*=-k*nu*L/Vr

Thus, since this is the driving force, one prescribe

Vr=k*nu*L

and Re= k*nu*L^2/D

But that assumes nu=constant, but if you have nu(u) the issue is different.

hunt_mat · August 10, 2023, 10:00

Quote:

Originally Posted by FMDenaro

In order to get an idea of the Reynolds number Vr*L/D, one needs to estimate Vr from the BCs:

du/dx=-k*nu -> du*/dx*=-k*nu*L/Vr

Thus, since this is the driving force, one prescribe

Vr=k*nu*L

and Re= k*nu*L^2/D

But that assumes nu=constant, but if you have nu(u) the issue is different.

Nu is separate to u, nu varies in space and time, it's the specific volume. Can you use the initial value of nu in the calculations?

FMDenaro · August 10, 2023, 11:13

Quote:

Originally Posted by hunt_mat

Nu is separate to u, nu varies in space and time, it's the specific volume. Can you use the initial value of nu in the calculations?

No, maybe the max value to get an idea of the max Re number

FMDenaro · August 10, 2023, 11:43

Quote:

Originally Posted by hunt_mat

Nu is separate to u, nu varies in space and time, it's the specific volume. Can you use the initial value of nu in the calculations?

Your M-files shows a solution not only oscillating but under instability, the values are aroud 10^73. Thus, a bug in the code or the numerical stability constraints are not satisfied.

FMDenaro · August 10, 2023, 11:48

A further big issue is that you have simply assigned

t=linspace(0,3,M);
x=linspace(0,1,N+1);
dt=t(2);
dx=x(2);

without any check about the CFDL and diffusive stability constraint

FMDenaro · August 10, 2023, 12:31

%Eulerian code for viscous Burgers equation with Neumann BC at x_max
clear;clc;
N=200; %This is the number of point in the x-variable, and also h.
M=1000; %This is the discretisation for the time variable;
t=linspace(0,0.0003,M);
x=linspace(0,1,N+1);
nu=0.1;
u=zeros(M,N+1); %This is the velocity;
dt=t(2);
dx=x(2);
D=0.01;
k=-0.1;
Re=1.0;

% initial condition for u
M=1
u0(1:N+1)=0.0 ;
%BC
u0(1)=0 ;
u0(N+1)=u0(N)-1 ;
u(M,1:N+1)=u0 ;

for k=1:1000,
% fluxes
for i=2:N+1,
conv=0.5*(u0(i)^2+u0(i-1)^2)/2.;
diff=(u0(i)-u0(i-1))/dx/Re;
f(i)=-conv+diff;
end

% Update
for i=2:N,
u1(i)=u0(i)+dt*(f(i+1)-f(i))/dx;
end
%BC
u1(1)=0;
u1(N+1)=u1(N)-1;

%Swap
u0(1:N+1)=u1(1:N+1);
u(k+1,1:N+1)=u1(1:N+1);

plot(u1)
end

pause

contour(u)
end

FMDenaro · August 10, 2023, 12:32

Se the few-lines code for solving the Burgers equation in FV formulation on a Eulerian grid.

Again, you have to compute the time step allowing the stability.

hunt_mat · August 11, 2023, 06:08

If I scale my variables according to the following:
$h=M\hat{h}$ , $t=\hat{t}/k$ , $\nu=\nu_{0}\hat{\nu}$ , and $u=Mk\nu_{0}\hat{u}$ leaves us with the following (cleaner) equations upon dropping the hats:
$\frac{\partial\nu}{\partial t}=\frac{\partial u}{\partial h}$
$\frac{\partial}{\partial t}(\nu u)=\frac{\partial }{\partial h}\left(\frac{\alpha}{\nu}\frac{\partial u}{\partial h}+\frac{u^{2}}{2}\right)$
The boundary conditions for u are

, $\frac{\partial u}{\partial h}\Bigg|_{h=1}=\nu(t,1)$ The constant $\alpha=\frac{D}{kM^{2}\nu_{0}^{2}}$
To do the stability analysis, do I look at the linear case?

FMDenaro · August 11, 2023, 06:23

Quote:

Originally Posted by hunt_mat

If I scale my variables according to the following:
$h=M\hat{h}$ , $t=\hat{t}/k$ , $\nu=\nu_{0}\hat{\nu}$ , and $u=Mk\nu_{0}\hat{u}$ leaves us with the following (cleaner) equations upon dropping the hats:
$\frac{\partial\nu}{\partial t}=\frac{\partial u}{\partial h}$
$\frac{\partial}{\partial t}(\nu u)=\frac{\partial }{\partial h}\left(\frac{\alpha}{\nu}\frac{\partial u}{\partial h}\right)$
The boundary conditions for u are

, $\frac{\partial u}{\partial h}\Bigg|_{h=1}=\nu(t,1)$ The constant $\alpha=\frac{D}{kM^{2}\nu_{0}^{2}}$
To do the stability analysis, do I look at the linear case?

The stability von Neumann analysis is always based on the linear equation.
In my code I used a very small time step as you can see. Try running.

hunt_mat · August 11, 2023, 07:52

Quote:

Originally Posted by FMDenaro

The stability von Neumann analysis is always based on the linear equation.
In my code I used a very small time step as you can see. Try running.

I like the fact you compute the fluxes as a separate function, I think that simplifies the code considerably, and it also allows the experimentation with different methods easily.

So I do a von Neumann stability analysis on the equations:
$\frac{\partial \nu}{\partial t}=\frac{\partial u}{\partial h}$
$\frac{\partial u}{\partial t}=\alpha\frac{\partial^{2}u}{\partial h^{2}}$
and find out the requirements for stability and try that?

FMDenaro · August 11, 2023, 08:42

Quote:

Originally Posted by hunt_mat

I like the fact you compute the fluxes as a separate function, I think that simplifies the code considerably, and it also allows the experimentation with different methods easily.

So I do a von Neumann stability analysis on the equations:
$\frac{\partial \nu}{\partial t}=\frac{\partial u}{\partial h}$
$\frac{\partial u}{\partial t}=\alpha\frac{\partial^{2}u}{\partial h^{2}}$
and find out the requirements for stability and try that?

A FVM is always coded by computing the fluxes sepatately.
For the stability analysis you have to work on the numerical schemes not on the linear PDEs.
Choose a scheme and the do that.

hunt_mat · August 11, 2023, 09:42

Quote:

Originally Posted by FMDenaro

A FVM is always coded by computing the fluxes sepatately.
For the stability analysis you have to work on the numerical schemes not on the linear PDEs.
Choose a scheme and the do that.

As I said, I'm relatively new to numerical stuff, and I only know selected bits, computing the fluxes separately has allowed my code to be easier to be read. I corrected my time step accordingly and it now converges beautifully.

Thank you for extending my knowledge of finite volume methods and also getting my code working.

hunt_mat · August 11, 2023, 10:06

So now I've got it working with the Euler method, are there any other methods I can use that will allow me to decrease the timestep?

FMDenaro · August 11, 2023, 11:35

Quote:

Originally Posted by hunt_mat

So now I've got it working with the Euler method, are there any other methods I can use that will allow me to decrease the timestep?

You have a lot of issues to learn about numerical methods. I strongly suggest to learn the fundamental topics in a good textbook.
For example, here you can read a lot of topics about numerical solution of PDEs

http://www.math.science.cmu.ac.th/do...20of%20PDE.pdf

hunt_mat · August 11, 2023, 16:42

Quote:

Originally Posted by FMDenaro

You have a lot of issues to learn about numerical methods. I strongly suggest to learn the fundamental topics in a good textbook.
For example, here you can read a lot of topics about numerical solution of PDEs

http://www.math.science.cmu.ac.th/do...20of%20PDE.pdf

I now have the thorny issue of replacing the constant D with a function of the porosity, and I have to think about how I relate the porosity to the specific volume.

FMDenaro · August 11, 2023, 16:55

Quote:

Originally Posted by hunt_mat

I now have the thorny issue of replacing the constant D with a function of the porosity, and I have to think about how I relate the porosity to the specific volume.

That is as same as happens in turbulence for eddy viscosity model. You have to compute D on any single face of the volumes by introducing the functional relation D=f(x,t).

hunt_mat · August 14, 2023, 07:16

A scaled version of the equations are:

$\frac{\partial\nu}{\partial t}=\frac{\partial u}{\partial h}$
$\frac{\partial u}{\partial t}=\frac{\partial}{\partial h}\left(\frac{\alpha\zeta(\theta)}{\nu}\frac{\partial u}{\partial h}\right)$
Where $\theta=1-\nu_{0}/\nu$ is the porosity, and the function is defined as
$\zeta(\theta)=\frac{2}{3\theta}(1-\theta)^{2}$
The boundary conditions are given as:
$u(t,0)=0,\quad\frac{\partial u}{\partial h}\Bigg|_{h=1}=-\nu(t,1)$
The code is playing up, and I was thinking it was something to do with the speed at which the diffusion happens, the speed seems to be increasing negatively quicker than the density seems to be increasing which I find odd.

Another thing I was thinking about is using the Neuman BC for u, to give me a BC for u at h=1, by writing the derivative as a one-sided stencil to get a value for u(t,1). Thoughts? Do you think that will things better?

FMDenaro · August 14, 2023, 10:06

Quote:

Originally Posted by hunt_mat

A scaled version of the equations are:

$\frac{\partial\nu}{\partial t}=\frac{\partial u}{\partial h}$
$\frac{\partial u}{\partial t}=\frac{\partial}{\partial h}\left(\frac{\alpha\zeta(\theta)}{\nu}\frac{\partial u}{\partial h}\right)$
Where $\theta=1-\nu_{0}/\nu$ is the porosity, and the function is defined as
$\zeta(\theta)=\frac{2}{3\theta}(1-\theta)^{2}$
The boundary conditions are given as:
$u(t,0)=0,\quad\frac{\partial u}{\partial h}\Bigg|_{h=1}=-\nu(t,1)$
The code is playing up, and I was thinking it was something to do with the speed at which the diffusion happens, the speed seems to be increasing negatively quicker than the density seems to be increasing which I find odd.

Another thing I was thinking about is using the Neuman BC for u, to give me a BC for u at h=1, by writing the derivative as a one-sided stencil to get a value for u(t,1). Thoughts? Do you think that will things better?

The Neumann Bc at x=xmax is what I used in the code I posted.

hunt_mat · August 14, 2023, 13:12

Quote:

Originally Posted by FMDenaro

The Neumann Bc at x=xmax is what I used in the code I posted.

I've been looking at your code, and I'm not sure what you've done. I don't quite understand how you've applied the BC. You've solved the viscous burgers equation, I am looking at a system, my conservation of mass, in terms of the finite volume approach is
$\int_{h_{j}-\frac{\delta h}{2}}^{h_{j}+\frac{\delta h}{2}}\frac{\partial \nu}{\partial t}dh=\left[ u\right]_{h_{j}-\frac{\delta h}{2}}^{h_{j}+\frac{\delta h}{2}}$
So on the end, I need to find an expression for $u_{h_{N}+\frac{\delta h}{2}}$
I think that If you use a one-sided backward derivative, you can use it to obtain it to get u on the right boundary.

August 14, 2023, 07:16	Final form of the equation	#78
hunt_mat Senior Member Matthew Join Date: Mar 2022 Location: United Kingdom Posts: 178 Rep Power: 4	A scaled version of the equations are: $\frac{\partial\nu}{\partial t}=\frac{\partial u}{\partial h}$ $\frac{\partial u}{\partial t}=\frac{\partial}{\partial h}\left(\frac{\alpha\zeta(\theta)}{\nu}\frac{\partial u}{\partial h}\right)$ Where $\theta=1-\nu_{0}/\nu$ is the porosity, and the function is defined as $\zeta(\theta)=\frac{2}{3\theta}(1-\theta)^{2}$ The boundary conditions are given as: $u(t,0)=0,\quad\frac{\partial u}{\partial h}\Bigg\|_{h=1}=-\nu(t,1)$ The code is playing up, and I was thinking it was something to do with the speed at which the diffusion happens, the speed seems to be increasing negatively quicker than the density seems to be increasing which I find odd. Another thing I was thinking about is using the Neuman BC for u, to give me a BC for u at h=1, by writing the derivative as a one-sided stencil to get a value for u(t,1). Thoughts? Do you think that will things better?

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
multiphaseEulerFoam FOAM FATAL IO ERROR	qutadah.r	OpenFOAM Running, Solving & CFD	11	December 10, 2021 20:18
SU2 7.0.7 Built on CentOS 7, parallel computation pyscript mpi exit error?	EternalSeekerX	SU2	3	October 9, 2020 18:28
Problem of simulating of small droplet with radius of 2mm	liguifan	OpenFOAM Running, Solving & CFD	5	June 3, 2014 02:53
Gradient evaluation in Finite Volume Methods	yidongxia	Main CFD Forum	7	August 6, 2012 10:23
Unstructured grid finite volume methods	Marcus	Main CFD Forum	3	December 5, 2000 00:25

August 10, 2023, 09:45		#62
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,831 Rep Power: 73	In order to get an idea of the Reynolds number VrL/D, one needs to estimate Vr from the BCs: du/dx=-knu -> du/dx=-knuL/Vr Thus, since this is the driving force, one prescribe Vr=knuL and Re= knuL^2/D But that assumes nu=constant, but if you have nu(u) the issue is different.

August 10, 2023, 11:48		#66
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,831 Rep Power: 73	A further big issue is that you have simply assigned t=linspace(0,3,M); x=linspace(0,1,N+1); dt=t(2); dx=x(2); without any check about the CFDL and diffusive stability constraint

August 10, 2023, 12:31		#67
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,831 Rep Power: 73	%Eulerian code for viscous Burgers equation with Neumann BC at x_max clear;clc; N=200; %This is the number of point in the x-variable, and also h. M=1000; %This is the discretisation for the time variable; t=linspace(0,0.0003,M); x=linspace(0,1,N+1); nu=0.1; u=zeros(M,N+1); %This is the velocity; dt=t(2); dx=x(2); D=0.01; k=-0.1; Re=1.0; % initial condition for u M=1 u0(1:N+1)=0.0 ; %BC u0(1)=0 ; u0(N+1)=u0(N)-1 ; u(M,1:N+1)=u0 ; for k=1:1000, % fluxes for i=2:N+1, conv=0.5(u0(i)^2+u0(i-1)^2)/2.; diff=(u0(i)-u0(i-1))/dx/Re; f(i)=-conv+diff; end % Update for i=2:N, u1(i)=u0(i)+dt(f(i+1)-f(i))/dx; end %BC u1(1)=0; u1(N+1)=u1(N)-1; %Swap u0(1:N+1)=u1(1:N+1); u(k+1,1:N+1)=u1(1:N+1); plot(u1) end pause contour(u) end

August 10, 2023, 12:32		#68
FMDenaro Senior Member Filippo Maria Denaro Join Date: Jul 2010 Posts: 6,831 Rep Power: 73	Se the few-lines code for solving the Burgers equation in FV formulation on a Eulerian grid. Again, you have to compute the time step allowing the stability.

August 11, 2023, 06:08		#69
hunt_mat Senior Member Matthew Join Date: Mar 2022 Location: United Kingdom Posts: 178 Rep Power: 4	If I scale my variables according to the following: $h=M\hat{h}$ , $t=\hat{t}/k$ , $\nu=\nu_{0}\hat{\nu}$ , and $u=Mk\nu_{0}\hat{u}$ leaves us with the following (cleaner) equations upon dropping the hats: $\frac{\partial\nu}{\partial t}=\frac{\partial u}{\partial h}$ $\frac{\partial}{\partial t}(\nu u)=\frac{\partial }{\partial h}\left(\frac{\alpha}{\nu}\frac{\partial u}{\partial h}+\frac{u^{2}}{2}\right)$ The boundary conditions for u are , $\frac{\partial u}{\partial h}\Bigg\|_{h=1}=\nu(t,1)$ The constant $\alpha=\frac{D}{kM^{2}\nu_{0}^{2}}$ To do the stability analysis, do I look at the linear case?

August 11, 2023, 10:06		#74
hunt_mat Senior Member Matthew Join Date: Mar 2022 Location: United Kingdom Posts: 178 Rep Power: 4	So now I've got it working with the Euler method, are there any other methods I can use that will allow me to decrease the timestep?