CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Sod Shock Tube solution using Nodal DGFEM 2D CNS Solver

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   July 30, 2016, 08:22
Angry Sod Shock Tube solution using Nodal DGFEM 2D CNS Solver
  #1
New Member
 
Quasar_89's Avatar
 
Suyash Sharma
Join Date: Feb 2016
Posts: 7
Rep Power: 10
Quasar_89 is on a distinguished road
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
Code:
rho  = 1;
rhou = y.^2;
rhov = 0;
Ener = (2*mu*x + pbar)/(gamma-1) + .5*(y.^4);


Sod's Shock Tube
Code:
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
%
Code:
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
Code:
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:

Code:
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
Attached Images
File Type: jpg 1D Sod's Solution.jpg (58.2 KB, 2 views)
File Type: png 2D_rho_sod.png (33.5 KB, 1 views)
File Type: png 2D_rhou_sod.png (52.0 KB, 1 views)
File Type: png 2D_Ener_sod.png (35.0 KB, 1 views)
File Type: png 2D_rho_sod_plot.png (14.7 KB, 2 views)

Last edited by Quasar_89; July 30, 2016 at 18:33.
Quasar_89 is offline   Reply With Quote

Reply

Tags
boundary condition, cns2d, dgfem, sod shock tube, time step size

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Shock tube simulation harish FLUENT 5 January 25, 2014 02:20
Working directory via command line Luiz CFX 4 March 6, 2011 20:02
Sod Shock Tube problem cfd Main CFD Forum 4 July 7, 2009 13:24
Need help in setting up a shock tube problem kalee FLUENT 0 November 2, 2007 14:11
Sod shock tube problem with real gas effects ganesh Main CFD Forum 0 October 17, 2006 01:22


All times are GMT -4. The time now is 15:12.