|
[Sponsors] |
![]() |
![]() |
#41 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,830
Rep Power: 73 ![]() ![]() ![]() |
Quote:
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. |
||
![]() |
![]() |
![]() |
![]() |
#42 | |
Senior Member
Matthew
Join Date: Mar 2022
Location: United Kingdom
Posts: 178
Rep Power: 4 ![]() |
Quote:
|
||
![]() |
![]() |
![]() |
![]() |
#43 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,830
Rep Power: 73 ![]() ![]() ![]() |
Quote:
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. |
||
![]() |
![]() |
![]() |
![]() |
#44 | |
Senior Member
Matthew
Join Date: Mar 2022
Location: United Kingdom
Posts: 178
Rep Power: 4 ![]() |
Quote:
What you say is probably okay for simple Eulerian co-ordinates. |
||
![]() |
![]() |
![]() |
![]() |
#45 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,830
Rep Power: 73 ![]() ![]() ![]() |
Quote:
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. |
||
![]() |
![]() |
![]() |
![]() |
#46 | |
Senior Member
Matthew
Join Date: Mar 2022
Location: United Kingdom
Posts: 178
Rep Power: 4 ![]() |
Quote:
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. |
||
![]() |
![]() |
![]() |
![]() |
#47 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,830
Rep Power: 73 ![]() ![]() ![]() |
Quote:
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? |
||
![]() |
![]() |
![]() |
![]() |
#48 | |
Senior Member
Matthew
Join Date: Mar 2022
Location: United Kingdom
Posts: 178
Rep Power: 4 ![]() |
Quote:
|
||
![]() |
![]() |
![]() |
![]() |
#49 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,830
Rep Power: 73 ![]() ![]() ![]() |
Quote:
|
||
![]() |
![]() |
![]() |
![]() |
#50 |
Senior Member
Matthew
Join Date: Mar 2022
Location: United Kingdom
Posts: 178
Rep Power: 4 ![]() |
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. |
|
![]() |
![]() |
![]() |
![]() |
#51 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,830
Rep Power: 73 ![]() ![]() ![]() |
Quote:
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. |
||
![]() |
![]() |
![]() |
![]() |
#52 | |
Senior Member
Matthew
Join Date: Mar 2022
Location: United Kingdom
Posts: 178
Rep Power: 4 ![]() |
Quote:
I think that was the big breakthrough I needed to understand in the code, as well as how to deal with the boundary conditions. |
||
![]() |
![]() |
![]() |
![]() |
#53 |
Senior Member
Lucky
Join Date: Apr 2011
Location: Orlando, FL USA
Posts: 5,728
Rep Power: 66 ![]() ![]() ![]() |
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. |
|
![]() |
![]() |
![]() |
![]() |
#54 |
Senior Member
Matthew
Join Date: Mar 2022
Location: United Kingdom
Posts: 178
Rep Power: 4 ![]() |
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:
Code:
%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); 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. |
|
![]() |
![]() |
![]() |
![]() |
#55 |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,830
Rep Power: 73 ![]() ![]() ![]() |
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. |
|
![]() |
![]() |
![]() |
![]() |
#56 | |
Senior Member
Matthew
Join Date: Mar 2022
Location: United Kingdom
Posts: 178
Rep Power: 4 ![]() |
Quote:
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. |
||
![]() |
![]() |
![]() |
![]() |
#57 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,830
Rep Power: 73 ![]() ![]() ![]() |
Quote:
Dirichlet condition at x=0 and non-homogeneous Neumann at x_max? What about k and D? |
||
![]() |
![]() |
![]() |
![]() |
#58 | |
Senior Member
Matthew
Join Date: Mar 2022
Location: United Kingdom
Posts: 178
Rep Power: 4 ![]() |
Quote:
I typed up the whole method and included it with this post. |
||
![]() |
![]() |
![]() |
![]() |
#59 | |
Senior Member
Filippo Maria Denaro
Join Date: Jul 2010
Posts: 6,830
Rep Power: 73 ![]() ![]() ![]() |
Quote:
That is relevant to get the stability constraints satisfied |
||
![]() |
![]() |
![]() |
![]() |
#60 |
Senior Member
Matthew
Join Date: Mar 2022
Location: United Kingdom
Posts: 178
Rep Power: 4 ![]() |
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.
|
|
![]() |
![]() |
![]() |
Tags |
finite volume method |
|
|
![]() |
||||
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 |