[Quasar_89]

Hello All,


I am working through my project involving DG-FEM code from Hesthaven & Warburton in their book on Nodal DG-FEM Alogorithms Analysis and Aplications.
http://github.com/tcew/nodal-dg/tree...Codes1.1/CFD2D


Related functions :

CurvedCNSdriver2D.m : Fintal Time and Polynomial Degree are edited in this driver code.

ChannelIC2D.m

ChannelBC2D.m

CurvedCNS2D.m
&
CurvedCNSdt2d.m

I’m trying to convert the case of channel flow into that of sod’s shock tube(closed ends) in CNS2D solver and for that I have tried to replace the ICs and the BCs with that of the Sod’s case. (Posted below : Channel Flow in Blue, Sod’s Case in Black).
http://www.csun.edu/~jb715473/exampl...1d.htm#density

I have kept otboxA01.neu associated with the channel flow as the grid. (13 nodes, 16 elements, triangulated on a 1x1 domain, all 4 boundaries as walls)
From what I have learnt, the case of no-slip boundary conditions shouldl have zero velocity at the boundaries at all times. Thus I’m trying to maintain the conserved momentum quantities at all 4 walls as that of the initial value zero, and the density and Energy stays constant on the 2 end walls.

With this setup the stable time step is not correctly computed and the dT value blows up pretty quickly for higher than 2nd order. This behavior is aggravated with higher polynomial order and vice-versa. I have been able to run the 1st and 2nd order cases by multiplying the dT value in CurvedCNSdT2D.m line-21 with .1 and .95 respectively. In fact in the 2nd order case the dT recovers after falling sharply.
The results with the below BC and IC for time = 0.2s are not evolving properly and in-fact the discontinuities are totally absent, and the results are pretty much under evolvoled.

I am not entirely sure where exactly the code is going wrong. I have not been able to find someone who has ran the code and known a fair bit about it. The problem here is limited time in my hands.




Initial Conditions : ChannelIC2D.m

Gamma = 1.4

Channel Flow

rho  = 1;
rhou = y.^2;
rhov = 0;
Ener = (2*mu*x + pbar)/(gamma-1) + .5*(y.^4);

Sod's Shock Tube

rho = ones(Np,K).*( (x<0.5) + 0.125*(x>=0.5));
rhou = zeros(Np,K);
rhov = zeros(Np,K);
Ener = ones(Np,K).*((x<0.5) + 0.1*(x>=0.5))/(gamma-1.0);

Boundary Conditions : ChannelBC2D.m

Channel Flow
%

rho(gauss.mapB)  = 1;
% rhou(gauss.mapB) = yB.^2;
% rhov(gauss.mapB) = 0;
% Ener(gauss.mapB) = (2*mu*xB + pbar)/(gamma-1) + .5*(yB.^4);

xBn = round(xB,4); % xB being the x co-ordinates of the gauss points on the boundary faces (mapB).

Sod's Shock Tube

rho(gauss.mapB(xBn==0)) = 1; rho(gauss.mapB(xBn==1))=0.125;
rhou(gauss.mapB) = 0;
rhov(gauss.mapB) = 0;
Ener(gauss.mapB(xBn==0)) = 1./(gamma-1); Ener(gauss.mapB(xBn==1))=0.1./(gamma-1);

Following is MATLAB code being used to plot the data contours:

xnew = 0:0.01:1;
ynew = xnew;
[Xnew,Ynew] = meshgrid(xnew,ynew);


for i =1:4
    aux = griddata(x,y,Q(:,:,i),Xnew(:),Ynew(:));
    Qnew(:,:,i) = reshape(aux,size(Xnew));
    figure
    contourf(Xnew,Ynew,Qnew(:,:,i),20,'LineColor','none')
    c = colorbar;
end

figure
contourf(Xnew,Ynew,sqrt(Qnew(:,:,2).^2 +Qnew(:,:,3).^2)./Qnew(:,:,1),20,'LineColor','none')
hold on
streamslice(Xnew,Ynew,Qnew(:,:,2)./Qnew(:,:,1),Qnew(:,:,3)./Qnew(:,:,1),10);

I need some advice on what exactly could be hampering the solution from evolving and also what exactly can lead to a stable time step size for higher orders. I also tried changing the RK stages in CurvedCNS2D.m