Skip to content

Commit 89d2214

Browse files
author
nianhuawong
committed
vanleer通量格式，检查，但是还是发散。
1 parent 83d7f71 commit 89d2214

File tree

1 file changed

+63
-41
lines changed

1 file changed

+63
-41
lines changed

‎src/Flux_Solver.cpp‎

Lines changed: 63 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -647,6 +647,7 @@ void Flux_Solver::VanLeer_Scheme()
647647

648648
void Flux_Solver::VanLeer_Scheme_X()
649649
{
650+
double nx = 1.0, ny = 0.0;
650651
VDouble fluxVector11(num_of_prim_vars);
651652
VDouble fluxVector22(num_of_prim_vars);
652653
for (int j = jst; j <= jed; j++)
@@ -659,9 +660,10 @@ void Flux_Solver::VanLeer_Scheme_X()
659660
double v1 = qField1[i][j][IV];
660661
double p1 = qField1[i][j][IP];
661662
double c1 = sqrt(fabs(gama * p1 / rho1));//声速取绝对值
662-
double V1 = sqrt(u1 * u1 + v1 * v1);
663-
double m1 = V1 / c1;
664-
663+
double vn1 = nx * u1 + ny * v1;
664+
double m1 = vn1 / c1;
665+
double h1 = Enthalpy(rho1, u1, v1, p1, gama);
666+
665667
double m1p, m2m;
666668
if (m1 >= 1)
667669
{
@@ -681,8 +683,9 @@ void Flux_Solver::VanLeer_Scheme_X()
681683
double v2 = qField2[i][j][IV];
682684
double p2 = qField2[i][j][IP];
683685
double c2 = sqrt(fabs(gama * p2 / rho2));//声速取绝对值
684-
double V2 = sqrt(u2 * u2 + v2 * v2);
685-
double m2 = V2 / c2;
686+
double vn2 = nx * u2 + ny * v2;
687+
double m2 = vn2 / c2;
688+
double h2 = Enthalpy(rho2, u2, v2, p2, gama);
686689

687690
if (m2 >= 1)
688691
{
@@ -701,37 +704,47 @@ void Flux_Solver::VanLeer_Scheme_X()
701704

702705
if (fabs(mn) < 1)
703706
{
704-
double fmass1 = rho1 * c1 * 0.25 * (m1 + 1) * (m1 + 1);
705-
double fmass2 = -rho2 * c2 * 0.25 * (m2 - 1) * (m2 - 1);
706-
707-
double term1 = pow((gama - 1) * V1 + 2.0 * c1, 2);
708-
double term2 = 0.5 * (u1 * u1 + v1 * v1 - V1 * V1);
709-
fluxVector11[IR] = fmass1;
710-
fluxVector11[IU] = fmass1 * ((-V1 + 2.0 * c1) / gama + u1);
711-
fluxVector11[IV] = fmass1 * ((-V1 + 2.0 * c1) / gama + v1);
712-
fluxVector11[IP] = fmass1 * ( term1 / 2.0 / (pow(gama, 2) - 1));
713-
714-
term1 = pow((gama - 1) * V2 - 2.0 * c2, 2);
715-
term2 = 0.5 * (u2 * u2 + v2 * v2 - V2 * V2);
716-
fluxVector22[IR] = fmass2;
717-
fluxVector22[IU] = fmass2 * ((-V2 - 2.0 * c2) / gama + u2);
718-
fluxVector22[IV] = fmass2 * ((-V2 - 2.0 * c2) / gama + v2);
719-
fluxVector22[IP] = fmass2 * (term1 / 2.0 / (pow(gama, 2) - 1));
707+
double fmass = rho1 * c1 * 0.25 * (m1 + 1) * (m1 + 1);
708+
double term0 = (-vn1 + 2.0 * c1) / gama;
709+
double term1 = pow((gama - 1) * vn1 + 2.0 * c1, 2);
710+
double term2 = 0.5 * (u1 * u1 + v1 * v1 - vn1 * vn1);
711+
//F+
712+
fluxVector11[IR] = fmass;
713+
fluxVector11[IU] = fmass * (u1 + nx * term0);
714+
fluxVector11[IV] = fmass * (v1 + ny * term0);
715+
//fluxVector11[IP] = fmass * h1;
716+
fluxVector11[IP] = fmass * (term1 / 2.0 / (pow(gama, 2) - 1) + term2);
717+
718+
//F-
719+
fmass = -rho2 * c2 * 0.25 * (m2 - 1) * (m2 - 1);
720+
term0 = (-vn2 - 2.0 * c2) / gama;
721+
term1 = pow((gama - 1) * vn2 - 2.0 * c2, 2);
722+
term2 = 0.5 * (u2 * u2 + v2 * v2 - vn2 * vn2);
723+
724+
fluxVector22[IR] = fmass;
725+
fluxVector22[IU] = fmass * (u2 + nx * term0);
726+
fluxVector22[IV] = fmass * (v2 + ny * term0);
727+
//fluxVector22[IP] = fmass * h2;
728+
fluxVector22[IP] = fmass * (term1 / 2.0 / (pow(gama, 2) - 1) + term2);
720729
}
721730
else if(mn >= 1)
722731
{
732+
//F+
723733
Inviscid_Flux_F(fluxVector11, rho1, u1, v1, p1);
734+
//F-
724735
fluxVector22[IR] = 0;
725736
fluxVector22[IU] = 0;
726737
fluxVector22[IV] = 0;
727738
fluxVector22[IP] = 0;
728739
}
729740
else if(mn <= -1)
730741
{
742+
//F+
731743
fluxVector11[IR] = 0;
732744
fluxVector11[IU] = 0;
733745
fluxVector11[IV] = 0;
734746
fluxVector11[IP] = 0;
747+
//F-
735748
Inviscid_Flux_F(fluxVector22, rho2, u2, v2, p2);
736749
}
737750

@@ -745,6 +758,7 @@ void Flux_Solver::VanLeer_Scheme_X()
745758

746759
void Flux_Solver::VanLeer_Scheme_Y()
747760
{
761+
double nx = 0.0, ny = 1.0;
748762
VDouble fluxVector11(num_of_prim_vars);
749763
VDouble fluxVector22(num_of_prim_vars);
750764
for (int i = ist; i <= ied; i++)
@@ -757,8 +771,9 @@ void Flux_Solver::VanLeer_Scheme_Y()
757771
double v1 = qField1[i][j][IV];
758772
double p1 = qField1[i][j][IP];
759773
double c1 = sqrt(fabs(gama * p1 / rho1));//声速取绝对值
760-
double V1 = sqrt(u1 * u1 + v1 * v1);
761-
double m1 = V1 / c1;
774+
double vn1 = nx * u1 + ny * v1;
775+
double m1 = vn1 / c1;
776+
double h1 = Enthalpy(rho1, u1, v1, p1, gama);
762777

763778
double m1p, m2m;
764779
if (m1 >= 1)
@@ -779,8 +794,9 @@ void Flux_Solver::VanLeer_Scheme_Y()
779794
double v2 = qField2[i][j][IV];
780795
double p2 = qField2[i][j][IP];
781796
double c2 = sqrt(fabs(gama * p2 / rho2));//声速取绝对值
782-
double V2 = sqrt(u2 * u2 + v2 * v2);
783-
double m2 = V2 / c2;
797+
double vn2 = nx * u2 + ny * v2;
798+
double m2 = vn2 / c2;
799+
double h2 = Enthalpy(rho2, u2, v2, p2, gama);
784800

785801
if (m2 >= 1)
786802
{
@@ -799,22 +815,28 @@ void Flux_Solver::VanLeer_Scheme_Y()
799815

800816
if (fabs(mn) < 1)
801817
{
802-
double fmass1 = rho1 * c1 * 0.25 * (m1 + 1) * (m1 + 1);
803-
double fmass2 = -rho2 * c2 * 0.25 * (m2 - 1) * (m2 - 1);
804-
805-
double term1 = pow((gama - 1) * V1 + 2.0 * c1, 2);
806-
double term2 = 0.5 * (u1 * u1 + v1 * v1 - V1 * V1);
807-
fluxVector11[IR] = fmass1;
808-
fluxVector11[IU] = fmass1 * ((-V1 + 2.0 * c1) / gama + u1);
809-
fluxVector11[IV] = fmass1 * ((-V1 + 2.0 * c1) / gama + v1);
810-
fluxVector11[IP] = fmass1 * (term1 / 2.0 / (pow(gama, 2) - 1));
811-
812-
term1 = pow((gama - 1) * V2 - 2.0 * c2, 2);
813-
term2 = 0.5 * (u2 * u2 + v2 * v2 - V2 * V2);
814-
fluxVector22[IR] = fmass2;
815-
fluxVector22[IU] = fmass2 * ((-V2 - 2.0 * c2) / gama + u2);
816-
fluxVector22[IV] = fmass2 * ((-V2 - 2.0 * c2) / gama + v2);
817-
fluxVector22[IP] = fmass2 * (term1 / 2.0 / (pow(gama, 2) - 1));
818+
double fmass = rho1 * c1 * 0.25 * (m1 + 1) * (m1 + 1);
819+
double term0 = (-vn1 + 2.0 * c1) / gama;
820+
double term1 = pow((gama - 1) * vn1 + 2.0 * c1, 2);
821+
double term2 = 0.5 * (u1 * u1 + v1 * v1 - vn1 * vn1);
822+
//G+
823+
fluxVector11[IR] = fmass;
824+
fluxVector11[IU] = fmass * (u1 + nx * term0);
825+
fluxVector11[IV] = fmass * (v1 + ny * term0);
826+
//fluxVector11[IP] = fmass * h1;
827+
fluxVector11[IP] = fmass * (term1 / 2.0 / (pow(gama, 2) - 1) + term2);
828+
829+
//G-
830+
fmass = -rho2 * c2 * 0.25 * (m2 - 1) * (m2 - 1);
831+
term0 = (-vn2 - 2.0 * c2) / gama;
832+
term1 = pow((gama - 1) * vn2 - 2.0 * c2, 2);
833+
term2 = 0.5 * (u2 * u2 + v2 * v2 - vn2 * vn2);
834+
835+
fluxVector22[IR] = fmass;
836+
fluxVector22[IU] = fmass * (u2 + nx * term0);
837+
fluxVector22[IV] = fmass * (v2 + ny * term0);
838+
//fluxVector22[IP] = fmass * h2;
839+
fluxVector22[IP] = fmass * (term1 / 2.0 / (pow(gama, 2) - 1) + term2);
818840
}
819841
else if (mn >= 1)
820842
{

0 commit comments

Comments
 (0)