ACSI Update: 2D Grids for Gas Dynamics Solver

Based on my extensive literature review, I have chosen a finite volume method with a structured grid for solving the equations of gas dynamics within a hybrid rocket motor. I will be using the Monotone Upstream-centered Scheme for Conservation Laws (MUSCL) in conjunction with AUSMPW+ approximate Riemann flux solver at the cell faces. These methods have demonstrated an enhanced capacity for capturing the oblique shocks and discontinuities present within a typical rocket motor and supersonic nozzle. Recently I have been investigating these methods in order to gain a more rigorous understanding of their function. You can download the documents I have read on my "Research Documents (CFD)" page under the heading "MUSCL Reconstructed AUSM Flux Splitting Scheme." So far, I have constructed a computational grid for an arbitrary radius function since the motors to be simulated are axisymmetric. The following diagram presents the theory behind the grid's construction. The code is also provided.

      PROGRAM NavierStokesSolver2D

      IMPLICIT NONE

      DOUBLE PRECISION Q(5, 4, 101, 101, 1001), n(4, 2, 101, 101),

     .L(4, 101, 101), A(101, 101)

     .X(5, 101), Y(5, 101, 101), T(1001)

      DOUBLE PRECISION Xmax, R(3, 101), Tmax, dx, dy(3, 101), dt

     .radius, height

      INTEGER i, j, k, Imax, Jmax, Kmax

     

      !Setting up the grids

      dt = Tmax / (Kmax-1)

      dx = Xmax / Imax

      DO i=1, Imax

          X(1, i) = (i - 1) * dx

          X(2, i) = (i - 1) * dx

          X(3, i) = (i - .5) * dx

          X(4, i) = i * dx

          X(5, i) = i * dx

         

          R(1, i) = radius(X(1, i))

          R(3, i) = radius(X(5, i))

          R(2, i) = (R(1, i) + R(3, i)) / 2

         

          dy(1, i) = 2 * R(1, i) / Jmax

          dy(2, i) = 2 * R(2, i) / Jmax

          dy(3, i) = 2 * R(3, i) / Jmax

         

          DO j=1, Jmax

              Y(1, i, j) = (j - 1) * dy(1, i) + height - R(1, i)

              Y(2, i, j) = j * dy(1, i) + height - R(1, i)

              Y(3, i, j) = (j - .5) * dy(2, i) + height - R(2, i)

              Y(4, i, j) = (j - 1) * dy(3, i) + height - R(3, i)

              Y(5, i, j) = j * dy(3, i) + height - R(3, i)

          END DO

      END DO

     

      !Finding cell areas, side lengths, and normal vectors

      DO i=1, Imax

          DO j=1, Jmax

          A(i, j) = dx * (dy(1, i) + dy(3, i)) / 2

          L(1, i, j) = dy(1, i)

          L(2, i, j) = sqrt((X(2,i) - X(3,i))**2 +

     .          (Y(2,i,j) - Y(3,i,j))**2)

          L(3, i, j) = dy(3, i)

          L(4, i, j) = sqrt((X(1,i) - X(4,i))**2 +

     .          (Y(1,i,j) - Y(4,i,j))**2)

             n(1, 1, i, j) = -1

             n(1, 2, i, j) = 0

             n(2, 1, i, j) = (Y(3,i,j) - Y(2,i,j)) / L(2, i, j)

             n(2, 2, i, j) = (X(2,i) - X(3,i)) / L(2, i, j)

             n(3, 1, i, j) = 1

             n(3, 2, i, j) = 0

             n(4, 1, i, j) = (Y(1,i,j) - Y(4,i,j)) / L(2, i, j)

             n(4, 2, i, j) = (X(4,i) - X(1,i)) / L(4, i, j)

          END DO

      END DO

         

      END PROGRAM NavierStokesSolver2D

   

      FUNCTION radius(x)

          IMPLICIT NONE

          DOUBLE PRECISION radius, x

          radius = 1

          RETURN

      END