TFG Chang + Gerard + Pau [Stokes, Navier-Stokes]

TFGChang/AeroForcesCalculation.m

@@ -0,0 +1 @@
+classdef AeroForcesCalculation < handle    properties (Access = public)        L        D        E    end    properties (Access = private)        mesh        nodesConditions        pressureFun    end    properties (Access = private)        connec        bMesh        pBoundary        normalVectors        lengthElement        centroid        centralPointsElement        nX        nY    end    methods (Access = public)        function obj = AeroForcesCalculation(cParams)            obj.init(cParams);            obj.identifyBoundaryEdges();            obj.createBoundaryMesh();            obj.assignPressureToBoundary();            obj.defineVariables();        end        function compute(obj)            obj.computeNormalVectors();            obj.convertNorVectorsToLanFunctions();            obj.computeDragForces();            obj.computeLiftForces();            obj.computeEfficiency();            %obj.plotResults();        end    end    methods (Access = private)        function init(obj,cParams)            obj.mesh              = cParams.mesh;            obj.nodesConditions   = cParams.nodesConditions;            obj.pressureFun       = cParams.pressureFun;        end        function identifyBoundaryEdges(obj)            obj.mesh.computeEdges();            e  = obj.mesh.edges.nodesInEdges;            bE = ismember(e,obj.nodesConditions);            bEIndices = find(prod(bE,2));            obj.connec = e(bEIndices,:);        end        function createBoundaryMesh(obj)            ss.coord    = obj.mesh.coord;            ss.connec   = obj.connec;            ss.kFace    = -1;            obj.bMesh   = Mesh.create(ss);            obj.bMesh   = obj.bMesh.computeCanonicalMesh();        end        function assignPressureToBoundary(obj)            obj.pBoundary = LagrangianFunction.create(obj.bMesh,1,obj.pressureFun.order);            pVals         = obj.pressureFun.fValues(obj.nodesConditions,1);            obj.pBoundary.setFValues(pVals);        end        function plotPressureDistribution(obj)            obj.pBoundary.plot();        end        function defineVariables(obj)            obj.normalVectors        = zeros(obj.bMesh.nelem,obj.bMesh.ndim);            obj.lengthElement        = zeros(obj.bMesh.nelem,1);            obj.centroid             = mean(obj.bMesh.coord);            obj.centralPointsElement = (obj.bMesh.coord(obj.bMesh.connec(:,1),:) + obj.bMesh.coord(obj.bMesh.connec(:,2),:))/2;        end        function computeNormalVectors(obj)            node1               = obj.bMesh.coord(obj.bMesh.connec(:,1),:);            node2               = obj.bMesh.coord(obj.bMesh.connec(:,2),:);            distanceNodes       = node2 - node1;            obj.lengthElement   = vecnorm(distanceNodes,2,2);            tanVector           = distanceNodes./obj.lengthElement;            obj.normalVectors   = tanVector * [0 1;-1 0];        end        function convertNorVectorsToLanFunctions(obj)            obj.nX = LagrangianFunction.create(obj.bMesh,1,'P0');            obj.nY = LagrangianFunction.create(obj.bMesh,1,'P0');            obj.nX.setFValues(obj.normalVectors(:,1));            obj.nY.setFValues(obj.normalVectors(:,2));        end        function computeDragForces(obj)            s.operation = @(x) -obj.pBoundary.evaluate(x).*obj.nX.evaluate(x);            s.mesh      = obj.mesh;            pNX         = DomainFunction(s);            obj.D       = Integrator.compute(pNX,obj.bMesh,2);        end        function computeLiftForces(obj)            s.operation = @(x) -obj.pBoundary.evaluate(x).*obj.nY.evaluate(x);            s.mesh      = obj.mesh;            pNY         = DomainFunction(s);            obj.L       = Integrator.compute(pNY,obj.bMesh,2);        end        function computeEfficiency(obj)            obj.E = obj.L/obj.D;        end        function plotResults(obj)            obj.plotPressureDistribution();            figure;            quiver(obj.centralPointsElement(:,1),obj.centralPointsElement(:,2),obj.normalVectors(:,1),obj.normalVectors(:,2));            hold on;            quiver(obj.centroid(1,1),obj.centroid(1,2),obj.D,0);            hold on;            quiver(obj.centroid(1,1),obj.centroid(1,2),0,obj.L);            hold on;            obj.bMesh.plot()            title("Lift and drag distribution");            xlabel("x");            ylabel("y");        end    endend

TFGChang/AeroForcesCalculationTestComputer.m

@@ -0,0 +1,78 @@
+classdef AeroForcesCalculationTestComputer < handle
+
+    properties (Access = public)
+        comparingLiftMessage
+        comparingDragMessage
+    end
+
+    properties (Access = private)
+        tolerance
+    end
+
+    properties (Access = private)
+        LiftOriginal
+        DragOriginal
+        LiftTest
+        DragTest
+    end
+
+    methods (Access = public)
+
+        function obj = AeroForcesCalculationTestComputer(cParams)
+            obj.init(cParams);
+        end
+
+        function compute(obj)
+            obj.compareLift();
+            obj.compareDrag();
+            result = [obj.comparingLiftMessage, newline, obj.comparingDragMessage];
+            disp(result);
+        end
+
+    end
+
+
+    methods (Access = private)
+
+        function init(obj,cParams)
+            obj.saveInput(cParams);
+            obj.loadOriginalfinalMesh();
+            obj.tolerance = 1e-4;
+        end
+
+        function saveInput(obj,cParams)
+            obj.LiftTest    = cParams.L;
+            obj.DragTest    = cParams.D;
+        end
+
+        function loadOriginalfinalMesh(obj)
+            load("E_GF_N2412.mat","Li","Di");
+            obj.LiftOriginal = Li;
+            obj.DragOriginal = Di;
+        end
+
+        function compareLift(obj)
+            try
+                relativeError = abs(obj.LiftOriginal - obj.LiftTest) / abs(obj.LiftOriginal);
+                assert(relativeError < obj.tolerance, ...
+                    sprintf('The lift force does not match the expected value.\nRelative error: %.4f', relativeError));
+                obj.comparingLiftMessage = 'The lift force corresponds to the expected value.';
+            catch ME
+                obj.comparingLiftMessage = ME.message;
+            end
+        end
+
+        function compareDrag(obj)
+            try
+                relativeError = abs(obj.DragOriginal - obj.DragTest) / abs(obj.DragOriginal);
+                assert(relativeError < obj.tolerance, ...
+                    sprintf('The drag force does not match the expected value.\nRelative error: %.4f', relativeError));
+                obj.comparingDragMessage = 'The drag force corresponds to the expected value.';
+            catch ME
+                obj.comparingDragMessage = ME.message;
+            end
+        end
+
+    end
+
+end

TFGChang/AirfoilOptimizer.m

@@ -0,0 +1,206 @@
+classdef AirfoilOptimizer < handle
+
+    properties (Access = public)
+        optimalParams
+    end
+
+    properties (Access = private)
+        features
+        learningRate
+        optimizer
+        tol
+        upperBC
+        lowerBC
+        ParamsMat
+        E
+        velFunMat
+        PFunMat
+    end
+
+    methods (Access = public)
+
+        function obj = AirfoilOptimizer(cParams)
+            obj.init(cParams);
+        end
+
+        function computeOptAirfoilParams(obj)
+            obj.computeOptimization();
+        end
+
+        function plotEEvolution(obj)
+            obj.plotE();
+        end
+
+        function generateAFSOPVideo(obj)
+            plotFun = @(s, i) AirfoilOptimizer.plotAirfoilContour(s, i);
+            VideoGenerator.compute("AirfoilOptimization", size(obj.ParamsMat,1),100, obj.ParamsMat, plotFun);
+        end
+
+        function generateVelVideo(obj)
+            plotFun = @(s, ~) TestNaca.plotVelocity(s);
+            VideoGenerator.compute('AirfoilOptimization-Velocity', size(obj.velFunMat,2), size(obj.velFunMat,2),obj.velFunMat, plotFun);
+        end
+
+        function generatePVideo(obj)
+            plotFun = @(s, ~) TestNaca.plotPressure(s);
+            VideoGenerator.compute('AirfoilOptimization-Pressure',  size(obj.PFunMat,2), size(obj.PFunMat,2),obj.PFunMat, plotFun);
+        end
+
+        function saveData(obj)
+            save("OptData.mat","obj");
+        end
+
+    end
+
+    methods (Access = private)
+
+        function init(obj,cParams)
+            obj.features                       = cParams.features;
+            obj.optimizer                      = cParams.optimizer;
+            obj.tol                            = cParams.tol;
+            obj.learningRate                   = cParams.learningRate;
+            obj.upperBC                        = cParams.upperBC;
+            obj.lowerBC                        = cParams.lowerBC;
+            obj.optimalParams                  = obj.features;
+            obj.E(1)                           = obj.optimizer.computeOutputValues(obj.optimalParams);
+            obj.ParamsMat(1,:)                 = obj.features;
+            [obj.velFunMat{1}, obj.PFunMat{1}] = computeVelPFun(obj);
+        end
+
+        function projected = projectParams(obj,params)
+            projected     = max(min(params, obj.upperBC), obj.lowerBC);
+            projected(2)  = projected(2) * (projected(1) > 1e-6);
+            projected(1)  = projected(1) * (projected(2) > 1e-6);
+        end
+
+        function computeOptimization(obj)
+            diff     = 1;
+            iter     = 1;
+            maxIter  = 101;%301
+
+            rho      = 0.9;   %0.9 %0.5
+            epsilon  = 1e-8;
+            cache    = zeros(size(obj.optimalParams));
+
+            while diff > obj.tol && iter < maxIter
+
+                % if (iter > 70)
+                %
+                %     obj.learningRate = 0.5;
+                % end
+                %
+                %  if (iter > 97)
+                %
+                %     obj.learningRate = 0.02;
+                %  end
+
+                gradient = obj.optimizer.computeGradient(obj.optimalParams);
+
+                cache = rho * cache + (1 - rho) * (gradient.^2);
+
+                adjustedGradient = obj.learningRate * gradient ./ (sqrt(cache) + epsilon);
+                updatedParams    = obj.optimalParams + adjustedGradient;
+
+                obj.optimalParams = obj.projectParams(updatedParams);
+
+                obj.E(end + 1) = obj.optimizer.computeOutputValues(obj.optimalParams);
+
+                diff = max(abs(obj.ParamsMat(end) - obj.optimalParams));
+
+                obj.ParamsMat(end + 1,:) = obj.optimalParams;
+
+                if mod(iter,2) == 0
+                    [obj.velFunMat{end + 1}, obj.PFunMat{end + 1}] = computeVelPFun(obj);
+                end
+
+                iter = iter + 1;
+            end
+        end
+
+        function plotE(obj)
+            figure;
+            plot(1:length(obj.E), obj.E);
+            xlabel('Iteration');
+            ylabel('Aerodynamic Efficiency');
+            title('Evolution of Airfoil Aerodynamic Efficiency vs Optimization Iteration');
+            grid on;
+        end
+
+        function [velFun,PFun] = computeVelPFun(obj)
+            Naca.flowType = "Stokes";
+            Naca.length   = 8;
+            Naca.height   = 4;
+            Naca.nx       = 420;
+            Naca.M        = obj.optimalParams(1);
+            Naca.p        = obj.optimalParams(2);
+            Naca.t        = obj.optimalParams(3);
+            Naca.chord    = 1;
+            Naca.AoA      = obj.optimalParams(4);
+
+            NacaClass = TestNaca(Naca);
+            NacaClass.compute();
+            velFun = NacaClass.velocityFun;
+            PFun   = NacaClass.pressureFun;
+
+        end
+
+
+    end
+
+    methods (Static)
+
+        function [xContour,yContour] = computeAirfoilContour(s)
+            x   = 0:0.001:1;
+            m   = s(1);
+            p   = s(2);
+            t   = s(3);
+
+            if (m > 1e-6)
+                yc   = (x<=p).*(m/p^2.*(2*p*x - x.^2)) + ...
+                    (x>p).*(m/(1-p)^2.*((1 - 2*p) + 2*p*x - x.^2));
+                dydx = (x<=p).*(2*m/p^2.*(p - x)) + ...
+                    (x>p).*(2*m/(1-p)^2.*(p - x));
+            else
+                yc   = 0;
+                dydx = 0;
+            end
+
+            yt   = 5*t.*(0.2969*sqrt(x)-0.1260*x-0.3516*x.^2+0.2843*x.^3-0.1036*x.^4);
+            theta = atan(dydx);
+
+            xu = x - yt.*sin(theta);
+            yu = yc + yt.*cos(theta);
+            xl = x + yt.*sin(theta);
+            yl = yc - yt.*cos(theta);
+
+            xContour = [xu, fliplr(xl)];
+            yContour = [yu, fliplr(yl)];
+        end
+
+        function [xRot,yRot] = rotateAirfoil(xContour,yContour,AoA)
+            xRot = xContour*cos(AoA) - yContour*sin(AoA);
+            yRot = xContour*sin(AoA) + yContour*cos(AoA);
+        end
+
+        function plotAirfoilContour(s,i)
+
+            AoA = -deg2rad(s(4));
+
+            [xContour,yContour] = AirfoilOptimizer.computeAirfoilContour(s);
+
+            [xRot,yRot] = AirfoilOptimizer.rotateAirfoil(xContour,yContour,AoA);
+
+            plot(xRot, yRot, 'k-', 'LineWidth', 1);
+            xlim([-0.2 1.2]);
+            ylim([-0.6 0.6]);
+            axis equal;
+            xlabel('x'); ylabel('y');
+            title(sprintf('Airfoil Shape Optimization - Iteration %d', i - 1));
+            grid on;
+
+        end
+
+
+    end
+
+end