When I call uq_display(myRBDO_RIA) I face the error below:
Error in uq_display_uq_rbdo (line 114)
xlim([0 - 0.05*(iter - 1), iter - 1 + 0.05*(iter - 1)])
Error in uq_module/Display (line 72)
this.displayFun(this);
Error in uq_display (line 11)
Display(module)
Error in uq_Rotor_Laval (line 65)
uq_display(myRBDO_RIA)
Please, could you support me in solving this issue?
Thank you,
Bárbara
Can you please provide a minimal reproducible example?
Best regards
Styfen
Hi Styfen,
Thank you very much for your support.
Bellow you find my example. Thank you.
File to run the RBDO:
clearvars
rng(100,'twister')
uqlab
%% COMPUTATIONAL MODEL
ModelOpts.mFile = 'Rotor_Descentrado_Mancal_Flexivel_reduced';
ModelOpts.isVectorized = false;
myModel = uq_createModel(ModelOpts);
%%% RELIABILITY-BASED DESIGN OPTIMIZATION (RBDO) SETUP
%% CONSTRAINT FUNCTION
%% HARD CONSTRAINT / LIMIT STATE FUNCTION: g_k(x)
RBDOOpts.LimitState.Model = myModel;
%% SOFT CONSTRAINT: f_j(x): Não se aplica.
%% PROBABILISTIC INPUT MODEL - ENVIRONMENTAL VARIABLES
Iopts.Marginals(1).Name = 'E' ;
Iopts.Marginals(1).Type = 'Gaussian';
Iopts.Marginals(1).Moments = [200e9 10e9];
myInput = uq_createInput(Iopts) ;
RBDOOpts.Input.EnvVar = myInput ;
%% DESIGN VARIABLES
RBDOOpts.Input.DesVar(1).Name = 'd';
% RBDOOpts.Input.DesVar(1).Type = 'constant';
RBDOpts.Input.DesVar(1).CoV = 0.05; % Coefficient of Variation - 5% of its mean or nominal value.
%% COST FUNCTION - Função a ser minimizada.
% c(d) : c é a função custo, e d são as varáveis de projeto
RBDOOpts.Cost.mString = '(pi * X(:,1).^2) ./ 4';
%% OPTIMIZATION SETUP
% DESIGN SPACE: the bounds of the design space are first defined as follows:
RBDOOpts.Optim.Bounds = [0.010; 0.016];
% Optionally, the starting point for the optimization algorithm can also be provided:
RBDOOpts.Optim.StartingPoint = 0.012;
% Define the target failure probability:
RBDOOpts.TargetPf = 0.01;
% Set the maximum number of iterations:
RBDOOpts.Optim.MaxIter = 500;
%% RELIABILITY-BASED DESIGN OPTIMIZATION (RBDO)
%% RIA - Reliability Index Approach
% RIA is a two-level approach with form in the inner loop.
% Use reliability index as constraint instead of the failure probability
RIAOpts = RBDOOpts;
RIAOpts.Type = 'RBDO';
RIAOpts.Method = 'RIA';
% RIAOpts.FORM.StartingPoint = [0.012];
% RIAOpts.Gradient.h = 0.001;
% Run the RBDO analysis:
myRBDO_RIA = uq_createAnalysis(RIAOpts);
% Print out a report of the reults:
uq_print(myRBDO_RIA)
% Display a graphical representation of the results:
uq_display(myRBDO_RIA)
File to run mfile: ‘Rotor_Descentrado_Mancal_Flexivel_reduced’:
function [gx] = Rotor_Descentrado_Mancal_Flexivel(X)
%% ENTRADA DE DADOS
%% SHAFT
length_shaft = 747e-3; % comprimento total do eixo [m]
%dia_shaft = 13e-3; % Diametro do eixo [m];
dia_shaft = X(:,1);
r_shaft = dia_shaft/2; % Raio do eixo [m]
area_shaft = pi * (dia_shaft/2).^2; % area do eixo [m^2]
I_shaft = (pi * dia_shaft.^4)/64; % Momento de inercia do eixo [kg*m^2]
% Conforme informado no slide da aula 8a da Prof Katia: Eixos de aco alpha ~ 0
beta = 1e-4; % Considerado beta = 1*10^(-4) conforme dado no slide da aula 7.
%% SHAFT MATERIAL: STEEL
E = X(:,2);
% E = 200e9; % Modulo de Elasticidade [Pa]
dens = 7850; % Densidade [kg/m^3]
Po = 0.3; % Coeficiente de Poisson do aco [N/A]
%% DISC
dia_disc = 90e-3; % diametro disco [m]
length_disc = 47e-3; % comprimento do disco [m]
mass_disc = dens * length_disc * ... % massa do disco [kg];
pi * (dia_disc/2)^2;
area_disc = pi * (dia_disc/2)^2; % area do disco [m^2]
volume_disc = area_disc * length_disc; % Volume do disco [m^3]
Ip_disc = (1 / 8) * mass_disc * (dia_disc^2); % Momento de inercia polar do disco [kg*m^2] formula conforme o livro do Friswell Tabela A1.1 do Apendice 1
Id_disc = Ip_disc / 2 + mass_disc * length_disc^2 / 12; % Momento de inercia diametral do disco [kg*m^2]
%% CILINDRICAL BEARING
h0 = 90e-6; % Folga radial do mancal (delta) [m]
length_brg = 20e-3; % Comprimento do mancal [m]
dia_brg = 30e-3; % Diametro do mancal [m]
r_brg = dia_brg/2; % Raio do mancal [m]
area_brg = (pi * dia_brg^2) / 4; % Area do mancal [m^2]
I_brg = (pi * dia_brg^4) / 64; % Momento de inercia da secao transversal do mancal [m^4]
brg_number = 2; % Numero de mancais
%% EXCITACAO
me = 3.7e-5; % Desbalanceamento de massa [Kg.m]
g = 9.81; % Aceleracao da gravidade [m/s^2]
%% OLEO
eta = 0.038233; %0.038233; % Viscosidade dinamica do Oleo ISO VG32 @ 30C [Pa.s]
%% ELEMENTOS FINITOS PARA 15 NÓS
% PARAMETROS DA MALHA
Nel = 14; % Numero de elementos;
Nnos = Nel + 1; % Numero de nos;
df = 4; % Numero de graus de liberdade
df_foundation = 2;
Nnos_foundation = 2;
a = 145e-3; % Distancia entre o mancal 1 (lado esquerdo) ate o centro do disco [m]
b = 435e-3; % Distancia entre o mancal 2 (lado direito) ate o centro do disco [m]
brg_span = a+b; % Distancia entre os centros dos mancais [m]
node_mancais = [3, 13]; % Indice dos nós dos mancais 3 e 15
node_disco = 6; % Indice do disco nó 6
mancais_ele = [2,3; 12,13]; % Elemento dos mancais A e B (para alterar as prop. devido ao diametro)
omegas = 1800;
n_omega = size(omegas, 2);
dists = [73.5, 10, 10, 67.5, 67.5, 70.833, 70.833, 70.833, 70.833, 70.833, 70.835, 10, 10, 73.5]*10^(-3); % Vetor que contem o comprimento no eixo X de cada elemento [m]
%% Vetores de Geometria do EIXO
area = zeros(Nel,1);
mom_ine = zeros(Nel,1);
dia = zeros(Nel,1);
area(:) = area_shaft; % Area do eixo
area(mancais_ele(1,:)) = area_brg; % Area do mancal 1
area(mancais_ele(2,:)) = area_brg; % Area do mancal 2
mom_ine(:) = I_shaft; % Mom. de Inercia do eixo
mom_ine(mancais_ele(1,:)) = I_brg; % Mom. de Inercia do mancal 1
mom_ine(mancais_ele(2,:)) = I_brg; % Mom. de Inercia do mancal 2
dia(:) = dia_shaft; % Diametro do eixo.
dia(mancais_ele(1,:)) = dia_brg; % Diametro do mancal 1.
dia(mancais_ele(2,:)) = dia_brg; % Diametro do mancal 2.
%% CALCULO DOS COEFICIENTES DOS MANCAIS
RA = mass_disc * g * b / brg_span; % Forca de reação do mancal A [N]
RB = mass_disc * g * a / brg_span; % Forca de reação do mancal B [N]
F0 = [RA; RB]; % Linha 1 = mancal A ; linha 2 = mancal B;
epsilon = zeros(length(node_mancais),length(omegas)); % Taxa de excentricidade [adimencional]
Fn = eta * length_brg^3 * omegas * r_brg / (2 * h0^2); % Força de referencia ou força de fricção
% Ss = F0 ./Fn; % Numero de Sommerfield Modificado (Conforme o livro do Kramer 6.30, pg 87
% Define initial guess
x0 = [h0, h0] ; % You can choose an initial guess here
%% Fsolve without Jacobian
options = optimoptions('fsolve', 'Display', 'off');
for j = 1:length(node_mancais)
for i = 1:length(omegas)
eps2 = eps^2;
S = @(eps) (pi/2) * (eps / ((1 - eps2)^2)) * sqrt(1 - eps2 + ((4 * eps / pi)^2)) - F0(j) / Fn(i);
inicial_guess = x0(j);
[epsilon_opt] = fsolve(S, inicial_guess, options);
epsilon(j,i) = epsilon_opt; % Store the result in the array
x0(j) = epsilon_opt; % Update the initial guess for the next iteration
end
end
%% COEFICIENTES PARA O MANCAL A E B
% Conforme Kramer 6.39 PDF pag.88
pi2 = pi^2;
epsilon2 = epsilon.^2;
epsilon4 = epsilon.^4;
A = 4 ./ (pi2 + (16 - pi2) .* epsilon2) .^ (3/2);
gama11 = (2 * pi2 + (16 - pi2) .* epsilon2) .* A;
gama12 = (pi / 4) .* ((pi2 - 2 * pi2 .* epsilon2 - (16 - pi2) .* epsilon4 ) ./ (epsilon .* sqrt(1 - epsilon2))) .* A;
gama21 = -(pi / 4) .* ((pi2 + (32 + pi2) .* epsilon2 + (32 - 2 * pi2) .* epsilon4) ./ (epsilon .* sqrt(1 - epsilon2))) .* A;
gama22 = ((pi2 + (32 + pi2) .* epsilon2 + (32 - 2 * pi2) .* epsilon4) ./ (1 - epsilon2)) .* A;
k11 = gama11 .* F0 ./ h0;
k12 = gama12 .* F0 ./ h0;
k21 = gama21 .* F0 ./ h0;
k22 = gama22 .* F0 ./ h0;
beta11 = (pi / 2) .* (sqrt(1 - epsilon2) ./ epsilon) .* (pi2 + (2 .* pi2 - 16) .* epsilon2) .* (A) ;
beta12 = -(2 * pi2 + (4*pi2 - 32) .* epsilon2 ) .* (A);
beta21 = beta12;
beta22 = (pi/2) .* ((pi2 + (48 - 2 * pi2) .* epsilon2 + (pi2 .* epsilon4)) ./ (epsilon .* (1 - epsilon2).^(1/2))) .* A;
d11 = beta11 .* F0 ./ (h0 .* omegas);
d12 = beta12 .* F0 ./ (h0 .* omegas);
d21 = beta21 .* F0 ./ (h0 .* omegas);
d22 = beta22 .* F0 ./ (h0 .* omegas);
%% Montagem da matriz global de todos os nos e elementos (sem condicao de contorno)
Kg = zeros(Nnos * df);
Mg = zeros(Nnos * df);
Gg = zeros(Nnos * df);
Cg = zeros(Nnos * df);
%% inci_eq: linha 1 = elemento 1: 4 graus de liberdade em cada no
ind_eq = zeros(Nel, 8);
for i=1:Nel
ind_eq(i,:) = 1 + ((i - 1) * 4) : (( i + 1 ) * 4);
end
for ii = 1:size(X(:,1))
for i = 1:Nel
he = dists(i);
% According to Kramer 16.31, pg 222
% Or according to Friswell 5.63, pg 189
%% Matriz de Massa do elemento
mi = dens * area(i);
c1 = 1 * 156 * mi * he / 420;
c2 = 1 * 54 * mi * he / 420;
c3 = he * 22 * mi * he / 420;
c4 = he * 13 * mi * he / 420;
c5 = he^2 * 4 * mi * he / 420;
c6 = he^2 * 3 * mi * he / 420;
me_upper = [c1, 0, 0, -c3, c2, 0, 0, c4;
0, c1, c3, 0, 0, c2, -c4, 0;
0, 0, c5, 0, 0, c4, -c6, 0;
0, 0, 0, c5, -c4, 0, 0, -c6;
0, 0, 0, 0, c1, 0, 0, c3;
0, 0, 0, 0, 0, c1, -c3, 0;
0, 0, 0, 0, 0, 0, c5, 0;
0, 0, 0, 0, 0, 0, 0, c5];
Me = me_upper + me_upper' - diag(diag(me_upper));
%% Matriz Giroscopica do elemento
mi_p = 2 * dens * mom_ine(i);
g1 = 6/5 * 1/he * mi_p;
g2 = 1/10 * 1 * mi_p;
g3 = 1/30 * he * mi_p;
g4 = 2/15 * he * mi_p;
ge_upper = [ 0, g1, g2, 0, 0, -g1, g2, 0;
0, 0, 0, g2, g1, 0, 0, g2;
0, 0, 0, g4, g2, 0, 0, g3;
0, 0, 0, 0, 0, g2, g3, 0;
0, 0, 0, 0, 0, g1, -g2, 0;
0, 0, 0, 0, 0, 0, 0, -g2;
0, 0, 0, 0, 0, 0, 0, g4;
0, 0, 0, 0, 0, 0, 0, 0];
Ge = ge_upper + ge_upper' - diag(diag(ge_upper));
%% Matriz de Rigidez do elemento
a1 = 12 * E * mom_ine(i) / he^3;
a2 = 6 * E * mom_ine(i) / he^2;
a3 = 4 * E * mom_ine(i) / he;
a4 = 2 * E * mom_ine(i) / he;
ke_upper = [ a1, 0, 0, -a2, -a1, 0, 0, -a2;
0, a1, a2, 0, 0, -a1, a2, 0;
0, 0, a3, 0, 0, -a2, a4, 0;
0, 0, 0, a3, a2, 0, 0, a4;
0, 0, 0, 0, a1, 0, 0, a2;
0, 0, 0, 0, 0, a1, -a2, 0;
0, 0, 0, 0, 0, 0, a3, 0;
0, 0, 0, 0, 0, 0, 0, a3];
Ke = ke_upper + ke_upper' - diag(diag(ke_upper));
%% Matriz de Amortecimento interno do elemento
Ce = beta * Ke;
%% Adicionando as matrizes dos elementos na matriz global
Kg(ind_eq(i,:),ind_eq(i,:)) = Kg(ind_eq(i,:), ind_eq(i,:)) + Ke;
Mg(ind_eq(i,:),ind_eq(i,:)) = Mg(ind_eq(i,:), ind_eq(i,:)) + Me;
Gg(ind_eq(i,:),ind_eq(i,:)) = Gg(ind_eq(i,:), ind_eq(i,:)) + Ge;
Cg(ind_eq(i,:),ind_eq(i,:)) = Cg(ind_eq(i,:), ind_eq(i,:)) + Ce;
end
end
%% ADICIONANDO OS COEFICIENTES DO MANCAL E A MASSA DO DISCO NOS RESPECTIVOS NÓS
Kgk = zeros(length(Kg),length(Kg),length(omegas));
Cgk = zeros(length(Cg),length(Cg),length(omegas));
Mgk = zeros(length(Mg),length(Mg),length(omegas));
Ggk = zeros(length(Gg),length(Gg),length(omegas));
n_var = size(Ggk, 1);
%% ADICIONANDO OS COMPONENTES DOS MANCAIS E DO DISCO:
for k = 1:n_omega
% Forcando copia
Kgk(:,:,k) = Kg(:,:);
Cgk(:,:,k) = Cg(:,:);
Mgk(:,:,k) = Mg(:,:);
Ggk(:,:,k) = Gg(:,:);
%% MANCAIS NA MATRIZ GLOBAL
for m = 1:brg_number
ind_brg = ind_eq(node_mancais(m), 1:4);
Kgk(ind_brg, ind_brg, k) = Kgk(ind_brg, ind_brg, k) + [k11(m, k), k12(m, k), 0, 0;
k21(m, k), k22(m, k), 0, 0;
0, 0, 0, 0;
0, 0, 0, 0];
Cgk(ind_brg, ind_brg, k) = Cgk(ind_brg, ind_brg, k) + [d11(m, k), d12(m, k), 0, 0;
d21(m, k), d22(m, k), 0, 0;
0, 0, 0, 0;
0, 0, 0, 0];
end
%% EFEITO DO DISCO NA MATRIZ GLOBAL:
ind_disc = ind_eq(node_disco, 1:4);
% Conforme o livro de Friswell (Equacao 3.49) ou Kramer 16.20, pg 219
Mgk(ind_disc, ind_disc, k) = Mgk(ind_disc, ind_disc, k) + diag([mass_disc, mass_disc, Id_disc, Id_disc]);
% Conforme o livro de Friswell (Equacao 3.49) ou Kramer 16.20, pg 219
Ggk(ind_disc, ind_disc, k) = Ggk(ind_disc, ind_disc, k) + [0, 0, 0, 0 ;
0, 0, 0, Ip_disc ;
0, 0, 0, 0 ;
0, -Ip_disc, 0, 0];
end
xi = zeros(Nnos * df, n_omega); % VETOR PARA GUARDAR O DESLOCAMENTO.
Fe_s = zeros(n_var, length(omegas)); % Para guardar o vetor na forma complexa.
for k = 1:n_omega
s = complex(0, omegas(k));
ind_disc = ind_eq(node_disco, 1:4);
Fe_s(ind_disc, k) = [ me * (s ^ 2), me * (s ^ 2), 0, 0];
Z(:,:,k) = (Mgk(:,:,k) .* s ^ 2) + ((Cgk(:,:,k) + (omegas(k) .* Ggk(:,:,k))) .* s) + Kgk(:,:,k);
xi(:, k) = Z(:,:,k) \ Fe_s(:,k);
end
glim = 0.00024;
gx = max(abs(xi(node_disco*df-3, :)));
gx = glim - gx;
end
Thanks for the nice reproducer! I can indeed reproduce the problem. The problem comes from the fact that the algorithm already converges after 1 iteration. UQLab then calculates the axis limits for the convergence curve (which is now a point) as [0, 0], which Matlab does not allow. This is something that UQLab could possibly handle better. If you need this plot, regardless of whether it makes sense, I can tell you how to patch it.
Hi @styfen.schaer ,
Thank you very much for your explanation, I will try to run this analysis considering a bigger range of omega (omegas = pi + (0:10*pi:2000)), thus I expect the analysis RIA will have more convergence points …
If I still face the problem I let you know.
Thank you a lot!
Best,