Steady vs unsteady solution code

src/Main.cpp

@@ -47,41 +47,64 @@ int main()
 	write_string_to_file("workdir/wing_nrm.dat", geom->normals_to_string());
 
 	PanelMethod panels;
-	double AoA = 5.0 * M_PI / 180.0;
-	double fvel = 5.0;
+	PanelMethod panels2;
+	double AoA = 1.0 * M_PI / 180.0;
+	double fvel = 0.001;
+
+	size_t NPANELS = 200;
+	double TSTEP = 0.01;
+
 	Eigen::Vector3d bvel = Eigen::Vector3d(0, fvel * std::sin(AoA), -fvel * std::cos(AoA));
 	Eigen::Vector3d omega = Eigen::Vector3d(0, 0, 0);
 
 	panels.thin_wings.push_back(geom);
 	panels.backward_difference_order = 1;
-	panels.shed_initial_wake(60, 0.01, bvel, omega);
+	panels.shed_initial_wake(NPANELS, TSTEP, bvel, omega);
 	panels.build_geometry_matrix();
 
+	panels2.thin_wings.push_back(geom);
+	panels2.shed_initial_wake(NPANELS, TSTEP, bvel, omega);
+	panels2.build_geometry_matrix();
+	panels2.build_dynamic(true);
+	panels2.solve(true);
+
 	// Initial solution that will be convected in the wake
 	panels.build_dynamic(true);
 	panels.solve(true);
 
 	panels.transfer_solution_to_wake();
 
 	double t = 0.0;
-	int MAX_IT = 60;
+	int MAX_IT = NPANELS;
 	Eigen::ArrayXd forces = Eigen::ArrayXd(MAX_IT);
+	Eigen::ArrayXd steady_forces = Eigen::ArrayXd(MAX_IT);
 	for(size_t i = 0; i < MAX_IT; i++)
 	{
-		AoA = 0.0 * M_PI / 180.0;
-		fvel = 5.0;
+		AoA = 1.0 * M_PI / 180.0;
+		fvel = 1.0;
 		bvel = Eigen::Vector3d(0, fvel * std::sin(AoA), -fvel * std::cos(AoA));
 		omega = Eigen::Vector3d(0, 0, 0);
 		panels.timestep(bvel, omega);
-		t += 1.0;
 		std::cout << "It: " << i + 1 << " / " << MAX_IT << std::endl;
 		panels.compute_cps(false);
-		forces(i) = panels.compute_aero_force()(1);
+		forces(i) = panels.compute_aero_force(true)(1);
+
+		panels2.shed_initial_wake(NPANELS, TSTEP, bvel, omega);
+		panels2.build_dynamic(true);
+		panels2.solve(true);
+		panels2.compute_cps(true);
+		steady_forces(i) = panels2.compute_aero_force(true)(1);
+		t += TSTEP;
+
 	}
 	std::stringstream s;
 	s << forces.transpose();
 	write_string_to_file("workdir/lift_over_time.dat", s.str());
 
+	std::stringstream s2;
+	s2 << steady_forces.transpose();
+	write_string_to_file("workdir/lift_over_time_steady.dat", s2.str());
+
 	//write_string_to_file("workdir/mat.dat", panels.geometry_matrix_to_string());
 	//write_string_to_file("workdir/dyn_mat.dat", panels.dynamic_matrix_to_string());
 	write_string_to_file("workdir/sln.dat", panels.solution_to_string(0));

src/PanelMethod.cpp

@@ -333,38 +333,39 @@ double PanelMethod::get_area(const ThinWing &wing, Eigen::Index panel)
 	return 0.5 * (v2 - v1).cross(v3 - v1).norm() + 0.5 * (v4 - v2).cross(v4 - v3).norm();
 }
 
-Eigen::Vector3d PanelMethod::compute_aero_force()
+Eigen::Vector3d PanelMethod::compute_aero_force(bool centerline)
 {
-	Vector3d acc; acc.setZero();
+	Vector3d acc;
+	acc.setZero();
 
-
-	/*for(size_t geom = 0; geom < thin_wings.size(); geom++)
+	if (centerline)
 	{
-		for(Index panel = 0; panel < thin_wings[geom]->quads.cols(); panel++)
+		for (Index i = 0; i < thin_wings[0]->num_chorwise - 1; i++)
 		{
-			Vector3d nrm = thin_wings[geom]->normals.col(panel);
-			double area = get_area(*thin_wings[geom], panel);
-
-			Vector3d force = nrm * cps(geom_sizes[geom] + panel);
-			acc += force * area;
+			Index idx = thin_wings[0]->num_spanwise / 2 * (thin_wings[0]->num_chorwise - 1) + i;
+			Vector3d nrm = thin_wings[0]->normals.col(idx);
+			double area = get_area(*thin_wings[0], idx);
+			Vector3d a = thin_wings[0]->vertices.col(thin_wings[0]->quads(0, idx));
+			Vector3d b = thin_wings[0]->vertices.col(thin_wings[0]->quads(3, idx));
+			double length = (b - a).norm();
+			Vector3d force = nrm * cps(idx);
+			acc += force * length;
 		}
-	}*/
-
-
-	// Center line lift
-
-	for(Index i = 0; i < thin_wings[0]->num_chorwise - 1; i++)
-	{
-		Index idx = thin_wings[0]->num_spanwise / 2 * (thin_wings[0]->num_chorwise - 1) + i;
-		Vector3d nrm =thin_wings[0]->normals.col(idx);
-		double area = get_area(*thin_wings[0], idx);
-		Vector3d a = thin_wings[0]->vertices.col(thin_wings[0]->quads(0, idx));
-		Vector3d b = thin_wings[0]->vertices.col(thin_wings[0]->quads(3, idx));
-		double length = (b - a).norm();
-		Vector3d force = nrm * cps(idx);
-		acc += force * length;
 	}
+	else
+	{
+		for(size_t geom = 0; geom < thin_wings.size(); geom++)
+		{
+			for(Index panel = 0; panel < thin_wings[geom]->quads.cols(); panel++)
+			{
+				Vector3d nrm = thin_wings[geom]->normals.col(panel);
+				double area = get_area(*thin_wings[geom], panel);
 
+				Vector3d force = nrm * cps(geom_sizes[geom] + panel);
+				acc += force * area;
+			}
+		}
+	}
 	return acc;
 }
 
@@ -582,8 +583,7 @@ void PanelMethod::compute_cps(bool STEADY)
 
 				// wake_scale is our timestep in essence
 				partial_phi /= std::pow(wake_scale, sln_hist.size() - 1);
-
-				cps(panel_idx) -= 2.0 * partial_phi;
+				cps(panel_idx) -= 2.0 * partial_phi / freestream.squaredNorm();
 			}
 		}
 	}

src/PanelMethod.h

@@ -106,7 +106,7 @@ class PanelMethod
 	void timestep(const Eigen::Vector3d& cur_vel, const Eigen::Vector3d& cur_angvel);
 
 	// Requires cps to be computed
-	Eigen::Vector3d compute_aero_force();
+	Eigen::Vector3d compute_aero_force(bool centerline);
 
 	std::string wake_geom_to_string(size_t for_geom);