[FMDenaro]

Originally Posted by hunt_mat

That's what I'm doing currently. I'm currently trying to figure out how to do the differencing for the mass differences.

Ok, let is start from the single Burgers equation.

In conservative form, the update is always given by

U(i,n+1) - U(i,n) = - dt * {F[u(i+1/2,n) -F[u(i-1/2,n)]}/h

Where U is the celle average of u. We disregard the difference and use directly u.

Now, introduce a vector F(i) with the attention of considering a staggered position of i-1/2.
First, you do a loop on F(i) for all the faces in the grid. Here is your specific method, upwind, centred and so on. After you have computed all the fluxes, the update is

u(i,n+1)=u(i,n) -(dt/h)*[F(i+1) - F(i)]

That’s all… of course, I assume you know all the required constraint fornthe numerical stability.

Try to use the first order upwind and the second order centred flux reconstruction starting from an initial sine velocity.

[hunt_mat]

Originally Posted by FMDenaro

Ok, let is start from the single Burgers equation.

In conservative form, the update is always given by

U(i,n+1) - U(i,n) = - dt * {F[u(i+1/2,n) -F[u(i-1/2,n)]}/h

Where U is the celle average of u. We disregard the difference and use directly u.

Now, introduce a vector F(i) with the attention of considering a staggered position of i-1/2.
First, you do a loop on F(i) for all the faces in the grid. Here is your specific method, upwind, centred and so on. After you have computed all the fluxes, the update is

u(i,n+1)=u(i,n) -(dt/h)*[F(i+1) - F(i)]

That’s all… of course, I assume you know all the required constraint fornthe numerical stability.

Try to use the first order upwind and the second order centred flux reconstruction starting from an initial sine velocity.

All this is very generic, I know this already. However, I don't think it's quite as simple as you seem to make it out. You have to explicitly write the fluxes in centred points of the cells, this takes some thought, especially when you have to think about derivatives with respect to h, as the mid points are tricky to get at.

[FMDenaro]

Originally Posted by hunt_mat

All this is very generic, I know this already. However, I don't think it's quite as simple as you seem to make it out. You have to explicitly write the fluxes in centred points of the cells, this takes some thought, especially when you have to think about derivatives with respect to h, as the mid points are tricky to get at.

As I wrote above, this is a very simple excercise that my under graduate students did in my course.
There is nothing complicated about the flux reconstruction, for example the upwind writes as u(i,n)^2/2 or u(i-1,n)^2/2, depending onnthe sign of the velocit. The diffusive flux is simply D*[u(i,n) - u(i-1,n)]/h. All is stored in F(i).
That’s all, no more than 10 lines of a matlab program.

[hunt_mat]

Originally Posted by FMDenaro

As I wrote above, this is a very simple excercise that my under graduate students did in my course.
There is nothing complicated about the flux reconstruction, for example the upwind writes as u(i,n)^2/2 or u(i-1,n)^2/2, depending onnthe sign of the velocit. The diffusive flux is simply D*[u(i,n) - u(i-1,n)]/h. All is stored in F(i).
That’s all, no more than 10 lines of a matlab program.

The difficulty is computing the values in the differences of dh at the mid points. Then the next issue is the boundary condition which you cannot simply insert.

What you say is probably okay for simple Eulerian co-ordinates.

[FMDenaro]

Originally Posted by hunt_mat

The difficulty is computing the values in the differences of dh at the mid points. Then the next issue is the boundary condition which you cannot simply insert.

What you say is probably okay for simple Eulerian co-ordinates.

Indeed I'm struggling in suggesting to quit from any idea of lagrangian approach and start solving directly your eqs.(1) in conservative form using a FVM. You can easily couple then each other.
Prescribing the BCs in terms of Dirichlet or Neumann conditions is not a real issue. Start first to set periodic condition, than if your code works fine you can insert different BCs.

[hunt_mat]

Originally Posted by FMDenaro

Indeed I'm struggling in suggesting to quit from any idea of lagrangian approach and start solving directly your eqs.(1) in conservative form using a FVM. You can easily couple then each other.
Prescribing the BCs in terms of Dirichlet or Neumann conditions is not a real issue. Start first to set periodic condition, than if your code works fine you can insert different BCs.

I understand what you're saying, I really do. However, there are two things that I am currently struggling with. What is the appropriate value I use for dh? My flux for one of my equations is (D/nu)*u_h+u^2/2. at h=0, I have u=0, how do I compute u_h at h=0: On the other boundary, I have u_h=k*nu, how do I compute u at h=1?

These are the issues I am currently struggling with. I can compute an updated h for the next timestep, and also the new co-ordinates by X_t=u.

[FMDenaro]

Originally Posted by hunt_mat

I understand what you're saying, I really do. However, there are two things that I am currently struggling with. What is the appropriate value I use for dh? My flux for one of my equations is (D/nu)*u_h+u^2/2. at h=0, I have u=0, how do I compute u_h at h=0: On the other boundary, I have u_h=k*nu, how do I compute u at h=1?

These are the issues I am currently struggling with. I can compute an updated h for the next timestep, and also the new co-ordinates by X_t=u.

Let me assume that the cell at i=1 (x=h/2) has the face i-1/2 at x=0 in your domain and requires the BC.
Your question assume that at x=0 is u_bc=0, right? Therefore the convective flux is zero, the diffusive flux at i-1/2 can be computed as



F(1) = D*[u(1)-u_bc]/(h/2)= 2*D*u(1)/h


at first order, or if you want a second order approximation, the diffusive flux can be computed by means of a quadratic lagrange interpolation between u_bc, u(1) and u(2).


What is the problem?

[hunt_mat]

Originally Posted by FMDenaro

Let me assume that the cell at i=1 (x=h/2) has the face i-1/2 at x=0 in your domain and requires the BC.
Your question assume that at x=0 is u_bc=0, right? Therefore the convective flux is zero, the diffusive flux at i-1/2 can be computed as



F(1) = D*[u(1)-u_bc]/(h/2)= 2*D*u(1)/h


at first order, or if you want a second order approximation, the diffusive flux can be computed by means of a quadratic lagrange interpolation between u_bc, u(1) and u(2).


What is the problem?

Just doing the derivative at the mid point? I still need to compute h there, would I do that as dh_i=int_{X_{i-1/2}}^{X_{i+1/2}}1/nudx, and then approximate the integral using the trapezium rule? Or do I note that from the conservation of mass, rho*dX=rho_0*dX_0, and so I can essentially use all my dh are independent of time?

[FMDenaro]

Originally Posted by hunt_mat

Just doing the derivative at the mid point? I still need to compute h there, would I do that as dh_i=int_{X_{i-1/2}}^{X_{i+1/2}}1/nudx, and then approximate the integral using the trapezium rule? Or do I note that from the conservation of mass, rho*dX=rho_0*dX_0, and so I can essentially use all my dh as constant?

Use the Eulerian constant step h

[hunt_mat]

Originally Posted by FMDenaro

Use the Eulerian constant step h

We're not doing Eulerian though. I think that if X_0 is the original co-ordinates at t=0, then I can write h as

h=\int_{0}^{X}(\rho(t,u)du=int_{0}^{X_0}rho_0(t,x) du

using the conservation of mass. If rho_0 is independent of space, ten dh will be constant. The physics of this is that the mass in each of the cells remains constant in time.

[FMDenaro]

Originally Posted by hunt_mat

We're not doing Eulerian though. I think that if X_0 is the original co-ordinates at t=0, then I can write h as

h=\int_{0}^{X}(\rho(t,u)du=int_{0}^{X_0}rho_0(t,x) du

using the conservation of mass. If rho_0 is independent of space, ten dh will be constant. The physics of this is that the mass in each of the cells remains constant in time.

Is your Eq.(1) valid? You have written the equation drho/dt= -d/dx (rho*u), thus why do you state that rho is constant? If you integrate in the Eulerian cell you see that the time variation of the cell-average value depends on the difference of the mass fluxes.
On the other hand, in lagrangian formulation, the density is not constant, not even along a path-line since Drho/Dt= -rho*du/dx.
What is always constant is based on the transport Reynolds theorem, that is the total mass in a material volume is constant in time.

This is the physics I see from Eqs.(1) you have written.

[hunt_mat]

Originally Posted by FMDenaro

Is your Eq.(1) valid? You have written the equation drho/dt= -d/dx (rho*u), thus why do you state that rho is constant? If you integrate in the Eulerian cell you see that the time variation of the cell-average value depends on the difference of the mass fluxes.
On the other hand, in lagrangian formulation, the density is not constant, not even along a path-line since Drho/Dt= -rho*du/dx.
What is always constant is based on the transport Reynolds theorem, that is the total mass in a material volume is constant in time.

This is the physics I see from Eqs.(1) you have written.

If the density is constant throughout the material at t=0, then dh will be constant. If however the density is specially variable, then dh will change from cell to cell, and you have to take account of that when you take the derivative u_h. It's not too hard to add this in though.

I think that was the big breakthrough I needed to understand in the code, as well as how to deal with the boundary conditions.

[LuckyTran]

You simply prescribe the BCs that are actually BCs and directly insert them. Whatever is left that you do not know is not a BC, it is a flux that needs to be determined the same as any other face flux that is not a boundary cell. Obviously you need to determine the face fluxes based on cell-centered values. But this is not a problem of how to apply a BC because it is not a BC. If u=a is your BC then likely u_h is not a BC but a variable you solve for. If u_h is a BC then likely u needs to be determined using the reconstruction.

Keep in mind that the face fluxs are not formally "the solution," only the values of the cell centers are the solution. Hence, you are free to determine the fluxes using any reconstruction scheme you desire, preferably accurate and stable ones, and still arrive at a solution. The problem is closed in cell values after you apply any reconstruction scheme. Just as an exercise, consider that a 0th order interpolation scheme for example will turn all face fluxes into their dependencies on the cell values. Obviously this is a terrible scheme in general (but actually it would work for diffusive fluxes) and you should at least use the cell gradients.

[hunt_mat]

So I coded up my system, and I still get the same errors, an oscillation in the velocity. This has been the main issue that been bugging me with previous approaches. The code I see used is:

%This code is meant to be the Lagrangian code for the sintering
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,3,M);
x=linspace(0,1,N+1);
nu=zeros(M,N); %This is the specific volume
u=zeros(M,N); %This is the velocity;
X=zeros(M,N+1); %This is the mass.
L=zeros(M,1);
d_var=zeros(1,N);
dt=t(2);
dx=x(2);
D=0.01;
k=-0.1;

%Give the initial values for the variables:
nu(1,:)=1;
rho_half=1./nu(1,:);
rho_0=zeros(1,N+1);
rho_0(2:N)=0.5*(rho_half(1,N-1)+rho_half(2:N));
rho_0(1)=1.5*rho_half(1)-0.5*rho_half(2);
rho_0(N+1)=1.5*rho_half(N)-0.5*rho_half(N-1);
h=cumtrapz(x,rho_0);
dh=diff(h);
X(1,:)=x;
L(1)=1;
%Do the time progression now:

for i=2:M
    nu(i,2:N-1)=nu(i-1,2:N-1)+dt*(u(i-1,3:N)-u(i-1,1:N-2))./(2*dh(2:N-1));
    A=(2*D./(nu(i-1,3:N)+nu(i-1,2:N-1))).*((u(i-1,3:N)-u(i-1,2:N-1))./(0.5*(dh(3:N)+dh(2:N-1))))+((u(i-1,3:N)+u(i-1,2:N-1)).^2)/8;
    B=(2*D./(nu(i-1,2:N-1)+nu(i-1,1:N-2))).*((u(i-1,2:N-1)-u(i-1,1:N-2))./(0.5*(dh(2:N-1)+dh(1:N-2))))+((u(i-1,2:N-1)+u(i-1,1:N-2)).^2)/8;
    d_var(2:N-1)=nu(i-1,2:N-1).*u(i-1,2:N-1)+dt./(dh(2:N-1)).*(A-B);
    %Add in the boundary conditions
    nu(i,1)=nu(i-1,1)+0.5*dt*(u(i-1,2)-u(i-1,1))/dh(1);
    nu(i,N)=nu(i-1,N)+dt*(u(i-1,N)-u(i-1,N-1))/dh(N);
    A=2*D/(nu(i-1,2)+nu(i-1,1))*(u(i-1,2)-u(i-1,1))/(0.5*(dh(1)+dh(2)))+((u(i-1,1)+u(i-1,2))^2)/8;
    B=-D*u(i-1,1)/(0.5*dh(1)*(1.5*nu(i-1,1)-0.5*nu(i-1,2)));
    d_var(1)=nu(i-1,1)*u(i-1,1)+(dt/dh(1))*(A-B);
    A=k*D+0.5*(1.5*u(i-1,N)-0.5*u(i-1,N-1));
    B=(2*D./(nu(i-1,N)+nu(i-1,N-1))).*((u(i-1,N)-u(i-1,N-1))./(0.5*(dh(N-1)+dh(N))))+((u(i-1,N-1)+u(i-1,N-1)).^2)/8;
    d_var(N)=nu(i-1,N)*u(i-1,N)+(dt/dh(N))*(A-B);
    u(i,:)=d_var./nu(i,:);
    %Compute the new grid
    X(i,2:N)=X(i-1,2:N)+0.5*dt*(u(i,1:N-1)+u(i,2:N));
    dum=dt*(1.5*u(i,N)-0.5*u(i,N-1));
    X(i,N+1)=X(i-1,N+1)+dt*(1.5*u(i,N)-0.5*u(i,N-1));
    L(i)=X(i,N+1);
end

plot(t,L);

When I have needed the values at the start and endpoints of the domain I have been using the approximation:
q_end=1.5*q_N-0.5*q_N-1 and a similar one for the initial point. I am at a loss at to what is going on.

[FMDenaro]

Please, address the initial and boundary conditions for u and rho. Also the lenght of the domain and the value of D.
I will try to solve that directly in the Eulerian framework.

[hunt_mat]

Originally Posted by FMDenaro

Please, address the initial and boundary conditions for u and rho. Also the length of the domain and the value of D.
I will try to solve that directly in the Eulerian framework.

The initial condition for u is 0 across the whole domain, and the initial condition for nu is that it's constant. The boundary conditions for u are u=0 at x=0, and u_h=-k*nu for X_max.
The length of the domain is initially equal to 1, the grid is defined by the equation dX/dt=u.

I've done the linear theory for this and I get the right type of results, so I'm at a loss.

[FMDenaro]

Originally Posted by hunt_mat

The initial condition for u is 0 across the whole domain, and the initial condition for nu is that it's constant. The boundary conditions for u are u=0 at x=0, and u_h=-k*nu for X_max.
The length of the domain is initially equal to 1, the grid is defined by the equation dX/dt=u.

I've done the linear theory for this and I get the right type of results, so I'm at a loss.

You missed to specify the initial and bcs for the system (1) you wrote.

Dirichlet condition at x=0 and non-homogeneous Neumann at x_max? What about k and D?

[hunt_mat]

Originally Posted by FMDenaro

You missed to specify the initial and bcs for the system (1) you wrote.

Dirichlet condition at x=0 and non-homogeneous Neumann at x_max? What about k and D?

Yes, correct about the boundary condition. D=0.01, k=0.1. I chose small values.

I typed up the whole method and included it with this post.

[FMDenaro]

Originally Posted by hunt_mat

Yes, correct about the boundary condition. D=0.01, k=0.1. I chose small values.

I typed up the whole method and included it with this post.

Later I will read, is better working in non dimensional form. However, the reference velocity must be deduced from the Neumann condition, what about nu ?

That is relevant to get the stability constraints satisfied

[hunt_mat]

Originally Posted by FMDenaro

Later I will read, is better working in non dimensional form. However, the reference velocity must be deduced from the Neumann condition, what about nu ?

That is relevant to get the stability constraints satisfied

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.