How to find the optimum initial value while with the convergence?

Hi, I have problem regarding with my system. I am trying to find the two optimum value of "ecci" and "mis_angle" that are within the range of some limits while converging the out put value values which are F_tot and M_tot to the F_app and M_app with the 10e-3 tolerance. However, I could not solve the problem or create looping system that satisfy the condition. I need to first converge the "load_tol=abs(F_app-F_tot)/F_tot; value" than, the mom_tol=abs(M_app-M_tot)/M_tot; however I could not manage the coding? Can someone help me to find the optimized value of ecci and mis_angle??

close all;
clc;
clear all;

%CALCULATION AND OPERATION PARAMETERS 
ni=79;                                      % number of node points in circumferential-direction 
nj=79;                                      % number of node points in L-direction 
hi=2*ni-1;                                  % number of height node’s in circumferential-direction
hj=2*nj-1;                                  % number of height node’s in L-direction
R=180e-3;                                   % raduis of the bearing (m)
Pout=0.0;                                   % pressure at outer boundary (top in z) (pa)
Pin=0.0;                                    % pressure at inner boundary (bottom in z)
Pini=0;                                     % inital guess at pressures (pa)
Pcav=0.0;                                   % cavitation pressure (pa)
C=0.14e-3;                                  % journal bearing clearance (m)
e_crit=1e-12;                               % convergence criteria
iteration=100000;                           % max number of iterations
ecci=0.8; %0.75:0.001:0.95;                                % eccentrcity ratio e/c
t=140e-3;                                   % 1/2 bearing length (m)
L=2*t;                                      % bearing length
mu=0.09;                                    % viscosity (Pa.s)
N_rpm=43;                                   % velocity (rpm)
omega=2*pi*N_rpm/60;                        % angular velocity (rad/sec)
U=omega*R;
B=L/R; 
F_app=2e5;                                  % applied load
M_app=4e3;                                  % applied moment

%MISASLIGNMENT PARAMETERS    
mis_angle=0.02; %0.01:0.001:0.03;           % misalignmend angle
ecci_star=tand(mis_angle)*L/C;              % eccentricity on the 
alfa=pi/2;
fi_0=pi/2;

%GRID PARAMETERS AND NON DIMENSIONAL PARAMETER                                 
% dimensional direction
Lx=2*pi*R; dx=Lx/(ni-1); x=0:dx/2:Lx;  % Circumferential direction
Ly=L; dy=Ly/(nj-1); y=0:dy/2:Ly;       % Length direction

% non-dimentional direction
Lxn=Lx/R; dxn=dx/R; xn=x/R;
Lyn=Ly/L; dyn=dy/L; yn=y/L;

%FILM THICKNESS OF FLUID MECHANISM  
[h_film] = film_misal(xn,yn,y,C,ecci,ecci_star,alfa,fi_0);

%FILM THICKNESS PARAMETERS FOR FURTHER CALCULATION
h_film_non=h_film./C;                     % non dimensionalized film thickness                                      

%PRESSURE SOLVER
[P_solve,P_solution,e_error]=Press_solv(ni,nj,Pini,dxn,dyn,h_film_non,mu,Pcav,B,omega,R,C,e_crit,iteration);

%PERFORMANCE PARAMETERS
%Load Carrying capacity, misalignment moment and attitude angle
[F_T_bar,F_R_bar,Fx,Fy,F_tot,att_angle,att_sig,act_angle,M_x,M_y,M_tot,M_ang]=Load_mom_cap(ni,nj,dxn,dyn,Lxn,Lyn,mu,omega,R,C,Ly,P_solve);
   
load_tol=abs(F_app-F_tot)/F_tot;
mom_tol=abs(M_app-M_tot)/M_tot;

function[h_film] = film_misal(xn,yn,y,C,ecci,ecci_star,alfa,fi_0)

h_film_func = @(xn,y)(1+ecci*cos(xn-fi_0)+ecci_star*(yn-0.5)'*cos(xn-alfa-fi_0));
h_film = h_film_func(xn,y)';
h_film=h_film*C;

end

function[P_solve,P_solution,e_error]=Press_solv(ni,nj,Pini,dxn,dyn,h_film_non,mu,Pcav,B,omega,R,C,e_crit,iteration)

for i=1:ni
    for j=1:nj
        P(i,j)=Pini;
    end
end

hhminus_i = zeros(size(P));
for i=2:ni-1
    for j=2:nj-1
        hhminus_i(i,j)= h_film_non(2*i-2,2*j-1);
    end
end
hhplus_i = zeros(size(P));
for i=2:ni-1
    for j=2:nj-1
        hhplus_i(i,j)= h_film_non(2*i,2*j-1);
    end
end
hhminus_j = zeros(size(P));
for i=2:ni-1
    for j=2:nj-1
        hhminus_j(i,j)= h_film_non(2*i-1,2*j-2);
    end
end
hhplus_j = zeros(size(P));
for i=2:ni-1
    for j=2:nj-1
        hhplus_j(i,j)= h_film_non(2*i-1,2*j);
    end
end

for i=1:1:ni
    for j=1:1:nj
        B_B(i,j)=(1/(dxn^2)*(hhminus_i(i,j)^3));
        A_A(i,j)=(1/(dxn^2)*(hhplus_i(i,j)^3));
        F_F(i,j)=(1/(dyn^2)*(hhminus_j(i,j)^3)*1/(B^2));
        E_E(i,j)=(1/(dyn^2)*(hhplus_j(i,j)^3)*1/(B^2));
        D_D(i,j)=B_B(i,j)+A_A(i,j)+F_F(i,j)+E_E(i,j);
        G_G(i,j)=((1/dxn)*(hhplus_i(i,j)-hhminus_i(i,j)));       
    end    
end
count=0;
for k=1:iteration
    count=count+1
    for ii=2:nj-1
        for jj=2:ni-1
            P(ii,jj)=1/D_D(ii,jj)*(B_B(ii,jj)*P(ii-1,jj)...
                                  +A_A(ii,jj)*P(ii+1,jj)...
                                  +F_F(ii,jj)*P(ii,jj-1)...
                                  +E_E(ii,jj)*P(ii,jj+1)...
                                  -G_G(ii,jj));              
            if P(ii,jj) <0
               P(ii,jj) =0;
            end
        end
    end
    e_error = (sum(P(:)) - sum(Pini(:)))/sum(P(:));
    Pini=P;
    if e_error< e_crit
        break;
    end

end
for ii=1:ni
    for jj=1:nj
        P_solve(ii,jj)=P(ii,jj);
    end
end

%Reynolds Cavitation Condition
for ii=1:ni
    for jj=1:nj
        if P_solve(ii,jj)<Pcav
           P_solve(ii,jj)=Pcav;
        end
    end
end

P_solution=P_solve.*6.*mu.*omega.*(R^2./C.^2);
end

function [F_T_bar,F_R_bar,Fx,Fy,F_tot,att_angle,att_sig,act_angle,M_x,M_y,M_tot,M_ang]=Load_mom_cap(ni,nj,dxn,dyn,Lxn,Lyn,mu,omega,R,C,Ly,P_solve)

x_phi=0:dxn:Lxn;
y_non=0:dyn:Lyn;
x_cos=cos(x_phi);
x_sin=sin(x_phi);

ni_even=ni-1; 
ni_odd=ni-2; 
nj_even=nj-1; 
nj_odd=nj-2;

for i=1:ni
    P_bar(i)=P_solve(i,1)+P_solve(i,nj);            % sum the first and last term
    for j=2:2:nj_even                               % add in the even terms
        P_bar(i)=P_bar(i)+2*P_solve(i,j);
    end
    for j=3:2:nj_odd                                % add in the odd terms
        P_bar(i)=P_bar(i)+4*P_solve(i,j);
    end
end

for i=1:ni
    R_bar=P_bar(1)*x_cos(1)+P_bar(ni)*x_cos(ni);	% sum the first and last term
    for i=2:2:ni_even                               % add in the even terms
        R_bar=R_bar+2*(P_bar(i)*x_cos(i));
    end
    for i=3:2:ni_odd                                % add in the odd terms
        R_bar=R_bar+4*(P_bar(i)*x_cos(i));
    end
end

for i=1:ni
    T_bar=P_bar(1)*x_sin(1)+P_bar(ni)*x_sin(ni);	% sum the first and last term
    for i=2:2:ni_even                     % add in the even terms
        T_bar=T_bar+2*(P_bar(i)*x_sin(i));
    end
    for i=3:2:ni_odd                      % add in the odd terms
        T_bar=T_bar+4*(P_bar(i)*x_sin(i));
    end
end

% Load carrying capacity with Simpson 1/3 Rule
F_R_bar=(dxn/3)*(dyn/3)*R_bar*6*omega*mu*(R^2/C^2)*Ly*R;
F_T_bar=(dxn/3)*(dyn/3)*T_bar*6*omega*mu*(R^2/C^2)*Ly*R;
F_tot=sqrt(F_R_bar^2+F_T_bar^2);

att_sig=atand(abs(F_T_bar/F_R_bar));
att_angle=(1-sign(F_R_bar))*90+sign(F_R_bar)*sign(F_T_bar)*att_sig;

%Acting angle
if F_R_bar>=0
   act_angle=180-asind(F_T_bar/F_tot);
end
if F_R_bar<=0
   act_angle=asind(F_T_bar/F_tot);
end

Fx=F_T_bar*cosd(act_angle)+F_R_bar*sind(act_angle);
Fy=F_T_bar*sind(act_angle)-F_R_bar*cosd(act_angle);
func_tht_ang=abs(Fx/Fy);

%%% Misalignment Moment
y_mom=y_non-0.5;
P_mom_solve=P_solve.*y_mom;
for i=1:ni
    M_bar(i)=P_mom_solve(i,1)+P_mom_solve(i,nj);            % sum the first and last term
    for j=2:2:nj_even                               % add in the even terms
        M_bar(i)=M_bar(i)+2*P_mom_solve(i,j);
    end
    for j=3:2:nj_odd                                % add in the odd terms
        M_bar(i)=M_bar(i)+4*P_mom_solve(i,j);
    end
end

for i=1:ni
    MT_bar=M_bar(1)*x_cos(1)+M_bar(ni)*x_cos(ni);	% sum the first and last term
    for i=2:2:ni_even                               % add in the even terms
        MT_bar=MT_bar+2*(M_bar(i)*x_cos(i));
    end
    for i=3:2:ni_odd                                % add in the odd terms
        MT_bar=MT_bar+4*(M_bar(i)*x_cos(i));
    end
end

for i=1:ni
    MR_bar=M_bar(1)*x_sin(1)+M_bar(ni)*x_sin(ni);	% sum the first and last term
    for i=2:2:ni_even                     % add in the even terms
        MR_bar=MR_bar+2*(M_bar(i)*x_sin(i));
    end
    for i=3:2:ni_odd                      % add in the odd terms
        MR_bar=MR_bar+4*(M_bar(i)*x_sin(i));
    end
end

M_x=(dxn/3)*(dyn/3)*MR_bar*6*omega*mu*(R^2/C^2)*Ly^2*R;
M_y=(dxn/3)*(dyn/3)*MT_bar*6*omega*mu*(R^2/C^2)*Ly^2*R;
M_tot=sqrt(M_x^2+M_y^2);

M_star=atand(abs(M_x/M_y));
M_ang=(1-sign(M_y))*90+sign(M_y)*sign(M_x)*M_star;
end

Answer 1

[Anurag Ojha]

Hello Sami

In order to find the optimized value of “ecci” and “mis_angle” you can use genetic algorithm. We can utilize the Genetic Algorithm to systematically explore the parameter space and identify the most effective values for "ecci" and "mis_angle" that minimize both load and moment errors.

The first step is to define an objective function that evaluates the performance of a system for a given set of parameters (“ecci” and “mis_angle”).

Adding a sample code for your reference:

function error = optimizationObjective(params, ni, nj, dxn, dyn, Lxn, Lyn, mu, omega, R, C, Ly, P_solve, F_app, M_app) 
ecci = params(1); 
mis_angle = params(2); 
ecci_star = tand(mis_angle) * L / C; 

% Call your existing functions to calculate F_tot and M_tot 
[~, ~, ~, ~, F_tot, ~, ~, ~, ~, ~, M_tot, ~] = Load_mom_cap(ni, nj, dxn, dyn, Lxn, Lyn, mu, omega, R, C, Ly, P_solve, ecci, ecci_star); 

% Calculate convergence errors 
load_tol = abs(F_app - F_tot) / F_tot; 
mom_tol = abs(M_app - M_tot) / M_tot; 

% Combine errors into a single objective value 
error = load_tol + mom_tol; 
end 

Now, use the Genetic Algorithm to search for optimal values of “ecci” and “mis_angle”.

% Define the optimization problem 
options = optimoptions('ga', 'Display', 'iter'); 
lb = [0.75, 0.01];  % Lower bounds for ecci and mis_angle 
ub = [0.95, 0.03];  % Upper bounds for ecci and mis_angle 

% Run the genetic algorithm 
[optimalParams, optimalError] = ga(@(params) optimizationObjective(params, ni, nj, dxn, dyn, Lxn, Lyn, mu, omega, R, C, Ly, P_solve, F_app, M_app), 2, [], [], [], [], lb, ub, [], options); 

% Display the optimal parameters 
disp('Optimal Parameters:'); 
disp(['ecci: ', num2str(optimalParams(1)), ', mis_angle: ', num2str(optimalParams(2))]); 
disp(['Optimal Error: ', num2str(optimalError)]); 

Adding MATLAB Documentation to get better understanding of genetic algorithm:

• https://www.mathworks.com/help/releases/R2021a/gads/ga.html?s_tid=doc_ta