PHI=atan2(V0,V2) ALPHA=THETA-PHI % calculate angles Incidentally, this leads us to our bemsolve function: Once we’ve done this and evaluated our cost function J and our gradient values, we return a single value for the cost function, and a single column gradient vector by repacking it. We need to unpack this vector so that we can use our a and b values to evaluate bemsolve. This is so that we can use MATLAB/Octave’s fminunc function to do the heavy lifting for us.
Write subscript gnu octave code#
Conveniently, we also need to do this in order to evaluate our cost function and our gradient, so rather than writing the same code out twice, we’ll just write another nested function to handle it.Ī second important point is that our a and b values are passed to bemfuncmin as a single vector. This our basic evaluation function that is, if we know a and b, we can just evaluate equations 1-5 and be done with it. J = sum() grad = % convert back to column vectorsĬlever readers will notice that this function is calling another function called bemsolve. =bemsolve(A0, B0) % solve system of equations using estimated A and B values Next, we have a nested function that calculates our cost function and gradient:Ī0 = ABs(1:size(BLADE,2)) B0 = ABs((size(BLADE,2)+1):end) % unpack the A0 and B0 vectorsĪ0 = A0' B0 = B0' % put into row vector form These values are all constants that we want to want to calculate outside of any loops we might need to implement. THETA = (pitch+BETA)*pi/180 % calculate theta along the bladeĪ0 = 0.1*ones(size(BLADE)) B0 = 0.01*ones(size(BLADE)) % initialize A0, B0 K3 = 4.0*pi*rho*v^2 K4 = K3*omega/v % calculate equation constants K1 = 0.5*rho*blades*CHORD K2 = 0.5*rho*blades*CHORD.*BLADE % calculate equation constants First, we calculate our constants and initialize memory for the operation:Įlement = BLADE(end)-BLADE(end-1) omega = (rpm/60)*2*pi % calculate element size and angular velocity So that’s the theory, let’s have a look at the implementation. We can then use these in conjunction with a more powerful solver like fmincg, or fminunc, which will let us find a solution for values that have difficulty converging. Where the subscript g denotes our initial guessed value, and subscript c denotes the iterated, calculated value. We construct two cost functions of the form:Īnd find the gradients dJ a/da, and dJ b/db: To get around this, I used a tactic from machine learning to leverage more powerful solution methods. This works after a fashion, but occasionally you run into values that won’t converge easily, even with relaxation factors.
We can then iterate by subbing in our new values until everything converges to within a satisfactory threshold.
The idea is that we initialize a guess for both a and b, calculate dT/dr, dQ/dr by rearranging equations 4 and 5, and then use our momentum theory equations (equations 8 and 9, rearranged in terms of a and b) to calculate our new a and b values. To get around this, we introduce Froude’s momentum theory, which gives us two equivalent expressions for dT and dQ that we can use to iterate for a and b: The difficulty in solving these equations is that we don’t know the axial induction factor a, and the swirl factor b, and we don’t have enough equations to solve for them. Given that we want to resolve the forces on our blade into normal and tangential components in the disc plane, we can show that:Īnd our C L and C D values are determined using the angle of attack alpha: We can see that for our propeller slice, our lift and drag functions are: Referring to the following diagram of a single blade element: Given that I have some experience working with and writing my own blade element solvers, I was roped into helping him out I also had the ulterior motive of wanting to improve my own understanding of the topic. The main goal of this was to look into an idea that my father had been mulling over, and to determine if it had any merit. So this is a quick write-up of a blade element solver I put together based on this code from an academic at Sydney University.
0 kommentar(er)
