Finite volume methods

FMDenaro · September 1, 2023, 15:05

Quote:

Originally Posted by hunt_mat

So yes, there is another equation that has to be solved, and yes, the equation for mass is indeed linear in this case, that's one of the advantages of this formulation. It's this equation however, that seems to be causing the problems when I'm solving the equations. The other equation can be written in the following way:
$u_{i+1,j}-u_{i,j}-\frac{\alpha\theta\delta t}{\delta h_{j}}\delta f_{i+1}-\frac{(1-\theta)\alpha\delta t}{\delta h_{j}}\delta f_{i}=0$

Where:
$\delta f=\frac{\zeta(\nu_{j+\frac{1}{2}})}{\nu_{j+\frac{1}{2}}}\frac{u_{i,j+1}-u_{i,j}}{(\delta h_{j+1}+\delta h_{j})/2}-\frac{\zeta(\nu_{j-\frac{1}{2}})}{\nu_{j-\frac{1}{2}}}\frac{u_{i,j}-u_{i,j-1}}{(\delta h_{j}+\delta h_{j-1})/2}$

The issue seems to be the Jacobian though.

I just had a further check, you could linearize the equation for u and solving for u(i+1) with a Thomas algorithm and then solve the mass equation where u(i+1) is known.

hunt_mat · September 1, 2023, 15:16

Quote:

Originally Posted by FMDenaro

Use a simple approach, this is my suggestion. The resulting accuracy is the same.

Isn't this the same as a one-iteration of Newton Raphson?

Is there a way of linearization numerically?

FMDenaro · September 1, 2023, 16:04

Quote:

Originally Posted by hunt_mat

Isn't this the same as a one-iteration of Newton Raphson?

Is there a way of linearization numerically?

I would simply use a first order Taylor expansion for ni

ni(i+1)=ni(i)+dt*dni/dt|i= ni(i)+dt*du/dx|i

hunt_mat · September 1, 2023, 16:49

Quote:

Originally Posted by FMDenaro

I would simply use a first order Taylor expansion for ni

ni(i+1)=ni(i)+dt*dni/dt|i= ni(i)+dt*du/dx|i

The problem is the derivative, in this case the Jacobian. The upper left hand corner is the identity matrix, and that's causing some issues. That's the problem here.

I will add that I tried to do this using ode15 and that crashed as well which was annoying.

FMDenaro · September 1, 2023, 16:53

Quote:

Originally Posted by hunt_mat

The problem is the derivative, in this case the Jacobian. The upper left hand corner is the identity matrix, and that's causing some issues. That's the problem here.

I will add that I tried to do this using ode15 and that crashed as well which was annoying.

The linearization would produce a simple linear algebric system, no Jacobian is required! In the matrix entries you have all functions that are known at time i. The matrix is tridiagonal and you solve with Thomas.

hunt_mat · September 1, 2023, 16:57

Quote:

Originally Posted by FMDenaro

The linearization would produce a simple linear algebric system, no Jacobian is required! In the matrix entries you have all functions that are known at time i. The matrix is tridiagonal and you solve with Thomas.

I don't care about Thomas, I use the / in matlab.

This is going to be a major change to the code. I would need to predefine the matrices wouldn't I? I can't simply define them numerically that I can using NR.

FMDenaro · September 1, 2023, 17:05

Quote:

Originally Posted by hunt_mat

I don't care about Thomas, I use the / in matlab.

This is going to be a major change to the code. I would need to predefine the matrices wouldn't I? I can't simply define them numerically that I can using NR.

It's a waste of computational resource and time. You have a tridiagonal matrix, memorized in only three vectors a,b,c. Stop. The Thomas algorithm is a two lines code. There is no sense in invoking the use of NR.

hunt_mat · September 4, 2023, 09:14

Quote:

Originally Posted by FMDenaro

If you substitute ui+1 in the previous one (mass), you get a single non-linear equation in the form

A(ni).nij+1 =known term

here you can introduce a linearization for A. Are you doing this way?

How would I go about this? Is there a quick way to obtain the matrix?

FMDenaro · September 4, 2023, 11:39

Quote:

Originally Posted by hunt_mat

How would I go about this? Is there a quick way to obtain the matrix?

The operator that one gets from the CN scheme applied on the diffusive term

I-(dt/Re)d^2/dx^2

can be inverted, for dt/Re, small as

[I-(dt/Re)d^2/dx^2]^-1=I+(dt/Re)d^2/dx^2 + ...

that would make explicit the solution.

hunt_mat · September 4, 2023, 12:33

Quote:

Originally Posted by FMDenaro

The operator that one gets from the CN scheme applied on the diffusive term

I-(dt/Re)d^2/dx^2

can be inverted, for dt/Re, small as

[I-(dt/Re)d^2/dx^2]^-1=I+(dt/Re)d^2/dx^2 + ...

that would make explicit the solution.

Only that's not the term is it? it's
$\frac{\partial}{\partial h}\left(\frac{\zeta(\nu)}{\nu}\frac{\partial u}{\partial h}\right)$
Are you suggesting that I don't use the FVM anymore with this?

LuckyTran · September 4, 2023, 13:03

Thomas algorithm is just the solver for the linear system as opposed to brute forcing it with gaussian elimination, SORm CLU, or using / in matlab. The linear system is the same whether you solve it with thomas or gaussian elimination or / or your wonky newton raphson. Ax=b doesn't suddenly not become Ax=b just because you change the solver.

If you are actually using / in matlab then you are done. But then you shouldn't also be doing Newton Raphson. Clearly you haven't built the proper linear system that should arise when you do FVM + implicit CN.

Trust me in the old days of matlab where it didn't support matrix division I also used to spam fsolve to solve linear systems because I was too lazy to code up my own linear solvers. It is a relic of the past. Let it be buried there.

FMDenaro · September 4, 2023, 13:48

Quote:

Originally Posted by hunt_mat

Only that's not the term is it? it's
$\frac{\partial}{\partial h}\left(\frac{\zeta(\nu)}{\nu}\frac{\partial u}{\partial h}\right)$
Are you suggesting that I don't use the FVM anymore with this?

I am focusing on the LHS of the equation, that is how you get the entries of the matrix. If you look at the LHS, you see something like

(I-alpha*A).u=q

for small alpha value, an approximation would be

u=(I+alpha*A+...).q

without solving the system you can get an expression for the solution.

But forget this method, just solve for u using the Thomas algorithm or any other solver you want.

hunt_mat · September 5, 2023, 06:44

Quote:

Originally Posted by FMDenaro

I am focusing on the LHS of the equation, that is how you get the entries of the matrix. If you look at the LHS, you see something like

(I-alpha*A).u=q

for small alpha value, an approximation would be

u=(I+alpha*A+...).q

without solving the system you can get an expression for the solution.

But forget this method, just solve for u using the Thomas algorithm or any other solver you want.

You've lost me now. What's your matrix A?

FMDenaro · September 5, 2023, 07:26

Quote:

Originally Posted by hunt_mat

Only that's not the term is it? it's
$\frac{\partial}{\partial h}\left(\frac{\zeta(\nu)}{\nu}\frac{\partial u}{\partial h}\right)$
Are you suggesting that I don't use the FVM anymore with this?

The matrix A is obtained when you discretize this term at time i+1

hunt_mat · September 12, 2023, 11:01

I seem to be having an error. To compute the new position of the interfacial points I solve:
$\frac{dX}{dt}=u$
But I get silly answers for the new values. Do I need to remap at each timestep? The remap involves the Jacobian, so how do I compute the new Jacobian? I'm now assuming that once I have computed the new Jacobian, I then solve:
$\frac{dX}{dX_{0}}=\frac{\rho_{0}}{\rho}$
to get the new proper co-ordinates?

FMDenaro · September 12, 2023, 11:32

Quote:

Originally Posted by hunt_mat

I seem to be having an error. To compute the new position of the interfacial points I solve:
$\frac{dX}{dt}=u$
But I get silly answers for the new values. Do I need to remap at each timestep? The remap involves the Jacobian, so how do I compute the new Jacobian? I'm now assuming that once I have computed the new Jacobian, I then solve:
$\frac{dX}{dX_{0}}=\frac{\rho_{0}}{\rho}$
to get the new proper co-ordinates?

I suppose your domain is shrinking, the shift of the interfacial point requires a remapping when the last cell reaches a negative volume. That could be easily, just modify only the cell where the interfacial point lies and let the others.

hunt_mat · September 12, 2023, 12:42

Quote:

Originally Posted by FMDenaro

I suppose your domain is shrinking, the shift of the interfacial point requires a remapping when the last cell reaches a negative volume. That could be easily, just modify only the cell where the interfacial point lies and let the others.

I fixed the problem. The key idea is that the mass in each cell is constant. I compute the density(in the form of the specific volume), and the new length of my cell is:
$dX=\nu dh$
If dx is the original length, I can compute the difference as dx-dX>0, and I simply take that away from the original length.

hunt_mat · September 14, 2023, 07:57

I've been running my code and I have a rather odd phenomenon in that the density. I have computed the fluxes at all interface points, including the boundary. The density at x=0 should remain the same for small times but it changes. There is no a priori no-flux condition at x=0, but it's not physical, and I don't understand what is going on exactly.

I have included a screen shot.

FMDenaro · September 14, 2023, 13:14

Quote:

Originally Posted by hunt_mat

I've been running my code and I have a rather odd phenomenon in that the density. I have computed the fluxes at all interface points, including the boundary. The density at x=0 should remain the same for small times but it changes. There is no a priori no-flux condition at x=0, but it's not physical, and I don't understand what is going on exactly.

I have included a screen shot.

Let me resume your problem to be sure.

1) At x=0 you have the condition u=0 but you do not fix a value for density?

2) At x=0 you have the face of the first finite volume, that is the first equations is resolved at x=dx/2, right?
3) Thus, the convective mass flux is zero at x=0 but mass will slowly increase because you have a negative velocity profile, so that mass enter from the right face at x=dx.
4) The lagrangian scheme you are using involves a diffusive flux at x=0?

hunt_mat · September 15, 2023, 04:22

Quote:

Originally Posted by FMDenaro

Let me resume your problem to be sure.

1) At x=0 you have the condition u=0 but you do not fix a value for density?

2) At x=0 you have the face of the first finite volume, that is the first equations is resolved at x=dx/2, right?
3) Thus, the convective mass flux is zero at x=0 but mass will slowly increase because you have a negative velocity profile, so that mass enter from the right face at x=dx.
4) The lagrangian scheme you are using involves a diffusive flux at x=0?

1) Correct. The flux for the conservation of mass involves u only.
2) Essentially, the first cell is at dh/2.
3) There is a diffusion equation for the velocity. The negative velocity is diffused in at h=1.
4) Yes.

September 12, 2023, 11:01	Computing the Eulerian position from the Lagrangian position	#135
hunt_mat Senior Member Matthew Join Date: Mar 2022 Location: United Kingdom Posts: 184 Rep Power: 4	I seem to be having an error. To compute the new position of the interfacial points I solve: $\frac{dX}{dt}=u$ But I get silly answers for the new values. Do I need to remap at each timestep? The remap involves the Jacobian, so how do I compute the new Jacobian? I'm now assuming that once I have computed the new Jacobian, I then solve: $\frac{dX}{dX_{0}}=\frac{\rho_{0}}{\rho}$ to get the new proper co-ordinates?

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

September 4, 2023, 13:03		#131
LuckyTran Senior Member Lucky Join Date: Apr 2011 Location: Orlando, FL USA Posts: 5,761 Rep Power: 66	Thomas algorithm is just the solver for the linear system as opposed to brute forcing it with gaussian elimination, SORm CLU, or using / in matlab. The linear system is the same whether you solve it with thomas or gaussian elimination or / or your wonky newton raphson. Ax=b doesn't suddenly not become Ax=b just because you change the solver. If you are actually using / in matlab then you are done. But then you shouldn't also be doing Newton Raphson. Clearly you haven't built the proper linear system that should arise when you do FVM + implicit CN. Trust me in the old days of matlab where it didn't support matrix division I also used to spam fsolve to solve linear systems because I was too lazy to code up my own linear solvers. It is a relic of the past. Let it be buried there.