Advanced Question Airfoils and rotation for capsule

asbjos

tuanibrO
Addon Developer
Joined
Jun 22, 2011
Messages
696
Reaction score
259
Points
78
Location
This place called "home".
Hello! I'm remaking vessel modules for the Mercury capsule, and having some trouble with the aerodynamic definitions.
I found coefficients for drag (cd), moment (cm) and lift (cl) in NTRS document 19710069960 (specifically pages 25-27).
I must admit: I'm not much used to aerodynamics, so my implementation could be wrong. But what I've done is to have the values for cd, cm and cl in the vertical lift function (vlift), and cd and cl in the horizontal lift (hlift), while setting cm = 0 in hlift, and halving the cd and cl from hlift and vlift, so that the sum of the cd is the original value from the NTRS document.

But, the problem is in the rotation. The cm orients the spacecraft blunt side first, as in real life. But historically, a 10 degree per second rotation was applied with roll thruster during reentry. During my simulation, the capsule simply refuses to roll. See video here
What changes do I need to do to make the capsule orient itself correctly when entering the atmosphere, but still allow rolling of the capsule?

The vlift and hlift functions are here:
PHP:
void ProjectMercury::vlift(VESSEL* v, double aoa, double M, double Re, void* context, double* cl, double* cm, double* cd)
{
	static const double cmp[13] = { // for mach 1
		0, -0.204, -0.193, -0.238, -0.337, -0.223, 0.03, 0.23, 0.31, 0.235, 0.194, 0.192, 0
	};

	static const double clp[13] = { // for mach 1
		0, -0.42, -0.04, 0.38, -0.38, -0.45, 0.06, 0.4, 0.24, -0.38, 0.02, 0.42, 0.0
	};

	// cd is dependent on mach, for aoa 0
	static const  double mach[14] = {
		0.0, 0.50, 0.7, 0.90, 1.00, 1.1, 1.30, 1.60, 2.00, 3.0, 5.00, 7.0, 9.6, 20.0
	};
	static const double cdp[14] = { // drag coeff at 0 AoA (blunt side first) for different mach numbers
		1.0, 1.02, 1.1, 1.23, 1.34, 1.4, 1.46, 1.49, 1.53, 1.6, 1.56, 1.5, 1.5, 1.5
	};

	double aoastep = 30.0 * RAD;
	aoa += PI;
	int idx = max(0, min(11, (int)(aoa / aoastep)));
	double d = aoa / aoastep - idx;

	int i = 0;
	while (i < 14 && M > mach[i])
	{
		i++;
	}

	double machFraction;
	if (i == 14)
	{
		*cd = cdp[13];
	}
	else if (i == 0)
	{
		*cd = cdp[0];
	}
	else
	{
		*cd = cdp[i - 1] + (cdp[i] - cdp[i - 1]) * (M - mach[i - 1]) / (mach[i] - mach[i - 1]);
	}

	*cl = clp[idx] + (clp[idx + 1] - clp[idx]) * d;
	*cm = cmp[idx] + (cmp[idx + 1] - cmp[idx]) * d;

	*cd *= 0.5;
}

void ProjectMercury::hlift(VESSEL* v, double aoa, double M, double Re, void* context, double* cl, double* cm, double* cd)
{
	static const double cmp[13] = {
		0, 0.204, 0.193, 0.238, 0.337, 0.223, -0.03, -0.23, -0.31, -0.235, -0.194, -0.192, 0
	};

	static const double clp[13] = {
		0, -0.42, -0.04, 0.38, -0.38, -0.45, 0.06, -0.45, -0.38, 0.38, -0.04, -0.42, 0.0
	};

	static const  double mach[14] = {
		0.0, 0.50, 0.7, 0.90, 1.00, 1.1, 1.30, 1.60, 2.00, 3.0, 5.00, 7.0, 9.6, 20.0
	};
	static const double cdp[14] = { // drag coeff at 0 AoA (blunt side first) for different mach numbers
		1.0, 1.02, 1.1, 1.23, 1.34, 1.4, 1.46, 1.49, 1.53, 1.6, 1.56, 1.5, 1.5, 1.5
	};

	int i = 0;
	while (i < 14 && M > mach[i])
	{
		i++;
	}

	double machFraction;
	if (i == 14)
	{
		*cd = cdp[13];
	}
	else if (i == 0)
	{
		*cd = cdp[0];
	}
	else
	{
		*cd = cdp[i - 1] + (cdp[i] - cdp[i - 1]) * (M - mach[i - 1]) / (mach[i] - mach[i - 1]);
	}

	double aoastep = 30.0 * RAD;
	aoa += PI;
	int idx = max(0, min(11, (int)(aoa / aoastep)));
	double d = aoa / aoastep - idx;
	*cl = clp[idx] + (clp[idx + 1] - clp[idx]) * d;
	*cm = 0.0;

	*cd *= 0.5;
}

The thrusters should have pressure-independent thrust/isp, and should give a 1 deg/s^2 acceleration. But also an initial roll is dampened by the airfoils.

So what have I done wrong?



Also, as a sidequestion: on page 27 of the NTRS document, at 0 deg angle and mach 0 the cd is listed as 1.00. Is it fair to assume that this is the actual cd, or do you imagine it has been normalised? I must say I haven't read the entire 111 page document, but haven't found a specific indication for normalisation so far.

---------- Post added 19th Nov 2019 at 08:49 ---------- Previous post was 18th Nov 2019 at 22:23 ----------

And this is how the airfoils are created, by the way:
PHP:
CreateAirfoil3(LIFT_VERTICAL, _V(0, 0, 0.635), vlift, NULL, 1.89, 1.89 * 1.89 * PI / 4.0, 1.0);
CreateAirfoil3(LIFT_HORIZONTAL, _V(0, 0, 0.635), hlift, NULL, 1.89, 1.89 * 1.89 * PI / 4.0, 1.0); // spherical symmetry
I.e. assuming capsule diameter 1.89 m and contact surface [math]\frac{\pi d^2}{4}[/math]

---------- Post added at 12:24 ---------- Previous post was at 08:49 ----------

I just discovered API call SetRotDrag, so I got the main question answered myself. :)
 

martins

Orbiter Founder
Orbiter Founder
Joined
Mar 31, 2008
Messages
2,448
Reaction score
462
Points
83
Website
orbit.medphys.ucl.ac.uk
If the capsule doesn't react to the thrust force to the extent expected, maybe the aerodynamic setup is too stiff? That is, the centre of pressure (the airfoil attack point) may be too far off the centre of gravity. In your code, the offset is 0.635. Is this value based on any data? What happens if you reduce it?
 

asbjos

tuanibrO
Addon Developer
Joined
Jun 22, 2011
Messages
696
Reaction score
259
Points
78
Location
This place called "home".
If the capsule doesn't react to the thrust force to the extent expected, maybe the aerodynamic setup is too stiff? That is, the centre of pressure (the airfoil attack point) may be too far off the centre of gravity. In your code, the offset is 0.635. Is this value based on any data? What happens if you reduce it?

Thank you, but I had failed to define a rotational drag, which caused this problem.
Now it rotates also under atmospheric flight, although I'm struggeling with real-life values for Project Mercury, or anything at all, really.

As for the offset, 0.635 was just a placeholder guesstimate. The actual historical value is 0.7239 metres.
 

kuddel

Donator
Donator
Joined
Apr 1, 2008
Messages
2,064
Reaction score
507
Points
113
Now it rotates also under atmospheric flight,
I have never compiled any aerodynamic model (vessel), so I might be wrong, but:
Shouldn't the rotations (lifts etc.) be direction depending? I mean in "forward" direction (launch) it should be minimal, while in backward (reentry) direction it should be maximal?!
 

asbjos

tuanibrO
Addon Developer
Joined
Jun 22, 2011
Messages
696
Reaction score
259
Points
78
Location
This place called "home".
I have never compiled any aerodynamic model (vessel), so I might be wrong, but:
Shouldn't the rotations (lifts etc.) be direction depending? I mean in "forward" direction (launch) it should be minimal, while in backward (reentry) direction it should be maximal?!

Good question. According to the documentation, the definition for the rot. drag coeffs. [math]r_{x,y,z}[/math] is [math]a_{x,y,z} = -\omega_{x,y,z} q S_y r_{x,y,z}[/math], with [math]a[/math] being angular deceleration, [math]\omega[/math] the angular velocity, [math]q[/math] the dynamic pressure, and [math]S_y[/math] the reference surface defined as the cross section projected along the y-axis.

While of course being a simplification, it does not hold for a dimensional analysis. If we assume the rotational drag coeffs to be dimensionless, we get the units [math]\frac{kg \cdot m \cdot rad}{s^3}[/math] for the angular deceleration.

Maybe MartinS could shed some light on this?
 
Last edited:

N_Molson

Addon Developer
Addon Developer
Donator
Joined
Mar 5, 2010
Messages
9,286
Reaction score
3,255
Points
203
Location
Toulouse
Same kind of issue here, I'm struggling simulating Tianwen's lander entry vehicle (a rather basic capsule) using a vertical and an horizontal airfoil. And that's a simple one, full ballistic, not supposed to generate any kind of lift.

Would be nice to have an example in the stock SDK package. The trick is that in the capsule case, the vessel isn't going forward. Its going backwards, as the heatshield is the "aft" part. Here's my stuff, my problem is that it spins madly on itself as soon as it gets some dynamic pressure :

C++:
void VLiftCoeff(VESSEL *v, double aoa, double M, double Re, void *context, double *cl, double *cm, double *cd)
{
    int i;
    const int nabsc = 9;
    static const double AOA[nabsc] = { -180 * RAD,-135 * RAD,-30 * RAD, 10 * RAD, 0 * RAD,10 * RAD,30 * RAD,135 * RAD,180 * RAD };
    static const double CL[nabsc] = { 1,      0.4,   0,      0,    0,     0,   0,     0.4,      1 };
    static const double CM[nabsc] = { 0,      0,  0.014, 0.0039, -0.006,-0.008,-0.010,     0,      0 };
    for (i = 0; i < nabsc - 1 && AOA[i + 1] < aoa; i++);
    double f = (aoa - AOA[i]) / (AOA[i + 1] - AOA[i]);
    *cl = CL[i] + (CL[i + 1] - CL[i]) * f;  // aoa-dependent lift coefficient
    *cm = CM[i] + (CM[i + 1] - CM[i]) * f;  // aoa-dependent moment coefficient
    double saoa = sin(aoa);
    double pd = 0.5 + 0.4*saoa*saoa;  // profile drag
    *cd = pd + oapiGetInducedDrag(*cl, 1.5, 0.7) + oapiGetWaveDrag(M, 0.75, 1.0, 1.1, 0.04);
    // profile drag + (lift-)induced drag + transonic/supersonic wave (compressibility) drag
}

// 2. horizontal lift component (vertical stabilisers and body)

void HLiftCoeff(VESSEL *v, double beta, double M, double Re, void *context, double *cl, double *cm, double *cd)
{
    int i;
    const int nabsc = 8;
    static const double BETA[nabsc] = { -180 * RAD,-135 * RAD,-90 * RAD,-45 * RAD,45 * RAD,90 * RAD,135 * RAD,180 * RAD };
    static const double CL[nabsc] = { 0,    +0.3,      0,   -0.3,  +0.3,     0,   -0.3,      0 };
    for (i = 0; i < nabsc - 1 && BETA[i + 1] < beta; i++);
    *cl = CL[i] + (CL[i + 1] - CL[i]) * (beta - BETA[i]) / (BETA[i + 1] - BETA[i]);
    *cm = 0.0;
    *cd = 0.5 + oapiGetInducedDrag(*cl, 1.5, 0.6) + oapiGetWaveDrag(M, 0.75, 1.0, 1.1, 0.04);
}

C++:
CreateAirfoil3(LIFT_VERTICAL, _V(0, 0, 0.1), VLiftCoeff, NULL, 1.7, 1.7 * 1.7 * PI, 1.0);
CreateAirfoil3(LIFT_HORIZONTAL, _V(0, 0, 0.1), HLiftCoeff, NULL, 1.7, 1.7 * 1.7 * PI, 1.0);
 
Last edited:
Top