Диссертация (1025147), страница 22

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 22)

Текст программы в среде Matlab для уточнения конечноэлементной модели, построенной в ПО PSE по результатам экспериментальногомодального анализаfunction Model_updatingclose all;clc;clear all;globalN_mat n_mat ro E_el mu n_mat_el n_freq_all...S_EN norm_S_EN freq_i dir_model dir_PSE N_clusters S_EN_0 N_el%%%%%%% SETTINGS %%%%%%%%%%% EXE and settings directory %%%%dir_PSE = 'C:\\Nikolaev\\Model_updating\\Model_updating_soft\\';%%%% Model directory %%%%dir_model'C:\\Nikolaev\\Model_updating\\Model_updating_soft\\PSE_model\\Cutter\\';%%%% Number of updating iterations %%%%N_iterations = 45;%%%% Maximum allowed dE in percents for whole proceduredE_max_global = 30;%%%% Maximum allowed dE in percents for one iterationdE_max = 15;%%%% Sensitivity ratio for multiple frequencies updatingS_EN_ratio_critical = 8;%%%% Number of parameters clustersN_clusters = 25;%%%% Number of converged iterations to stop the solutionN_iterations_converged = 3;%%%% Minimal error in percents for convergence check %%%%error_critical = 2;%%%% Number of updated frequencies %%%%N_freq_update = 1;%%%% Reading test results %%%%freq_ref = load([dir_model '\test_freq.txt']); freq_ref = freq_ref';n_freq_all = length(freq_ref);%%%% Check the number of updated frequencies %%%%if (N_freq_update > length(freq_ref))N_freq_update = length(freq_ref);endk_convergence = 0;%%%%% Updating procedure %%%%for it_i = 1:N_iterations=175%%%% Reading FE model, PSE launch %%%%f = fopen([dir_model '\\ans_model_' num2str(it_i-1) '.txt'], 'r');PSE_model_import(f, it_i)%%%% Reading frequencies and sensitivity matrix %%%f_results = fopen([dir_model '\EIGF_energy.bin'], 'r');PSE_results_read(f_results, it_i)%%% Reference and currecnt frequencies for monitoring %%%freq_i_all = freq_i;freq_ref_all = freq_ref;%%% Reference and currecnt frequencies for updating %%%freq_i_update = freq_i(1:N_freq_update);freq_ref_update = freq_ref(1:N_freq_update);%%%% Initial error, initial E %%%%if (it_i == 1)E0 = E_el(1);err_freq_0 = abs((freq_i_update - freq_ref_update)./freq_ref_update);endfreq_i%%%% Current frequency errorerr_freq_i = 100*abs((freq_i_update - freq_ref_update)./freq_ref_update)freq_error(:,it_i) = 100*((freq_i_all - freq_ref_all)./freq_ref_all);freq_FE(:,it_i) = freq_i_update;%%%% Convergence check %%%%length(find(err_freq_i < error_critical))if (length(find(err_freq_i < error_critical)) == N_freq_update)k_convergence = k_convergence + 1;disp('Iteration has been converged')endif (k_convergence == N_iterations_converged)it_i = N_iterations;disp('Solution has been converged')return;end% %%%% Sensitivity matrix regularization %%%%S_EN = S_EN_0(:,1:N_freq_update);%%% Sensitivity normalization %%%%% for i = 1:N_freq_update%norm_S_EN(i,1) = norm(S_EN_0(:,i),2);%S_EN(:,i) = S_EN_0(:,i)./norm_S_EN(i,1);% end%%%% Sensitivity separation %%%%for i = 1:length(S_EN(:,1))if (S_EN(i,1)/S_EN(i,2) < S_EN_ratio_critical)S_EN(i,:) = 0;end[S_EN_max_i, ind_max_i] = max(S_EN(i,:));S_EN(i,ind_max_i) = 0;176S_EN_submax_i = max(S_EN(i,:));S_EN_ratio = S_EN_max_i/S_EN_submax_i;S_EN(i,:) = 0;if (S_EN_ratio > S_EN_ratio_critical)S_EN(i,ind_max_i) = S_EN_max_i*norm_S_EN(ind_max_i);endendfor i = 1:N_freq_updateif (err_freq_i(i) < error_critical)disp([num2str(i) ' frequency is ok'])endend%%%% Least square solution %%%%T_pinv = pinv(S_EN);%%%% dE calculation %%%d_freq_squared = freq_ref_update.^2 - freq_i_update.^2;dE = T_pinv'*d_freq_squared;%%%% Initial dE %%%%if (it_i == 1)dE_max_it_1 = max(max(abs(dE)));end%%%% Scalingscale_factor = (dE_max/100)*(E0/dE_max_it_1);dE = dE*scale_factor;%%%% dE constraintsE_min = E0*(1 - 0.01*dE_max_global);E_max = E0*(1 + 0.01*dE_max_global);%%%%%% Just negative dE are permitted %%dE(dE>0) = -1*dE(dE>0);%%%% Parameters clustering %%%%[ind_cluster, dE_cluster] = kmeans(dE,N_clusters);for i = 1:N_clustersind_2 = ind_cluster==i;dE(ind_2) = dE_cluster(i,1);endif (N_mat~=N_clusters)E_el = E_el(1)*ones(N_clusters,1);E_el = E_el';endn_mat = 1:N_clusters;ro = ro(1)*ones(N_clusters,1);mu = mu(1)*ones(N_clusters,1);%%%% Parameters updating %%%%n_mat_el = ind_cluster;E_el = E_el + dE_cluster';177E_el(E_el < E_min) = E_min;E_el(E_el > E_max) = E_max;%%%% Model export %%%%f = fopen([dir_model '\\ans_model_' num2str(it_i) '.txt'], 'w');PSE_model_export(f);figure(2)hold on; grid on;plot(freq_i_update, 'bo-')plot(freq_ref_update, 'ro-')figure(3)hold on; grid on;for i = 1:n_freq_allplot(freq_error(i,:), 'bo-')endsave([dir_model 'Convergence.mat'], 'freq_FE', 'freq_ref');end%%%% Convergence trackingfigure(3)hold on; grid on;plot(freq_error, 'b-')end%%%% Reading FE model, PSE launch %%%%function PSE_model_import(f, it_i)global n_DOF N_nodes N_el scale_u scale_F N_mat...n_mat ro E_el mu N_el_type el_type rot_vector node_num...cord n_type_el n_mat_el some_num el_number ind_el n_fixed_nodes...fixations n_freq_all dir_PSE dir_modeln_DOF = fscanf(f,'%d',1);N_nodes = fscanf(f,'%d',1);N_el = fscanf(f,'%d',1);scale_u = fscanf(f,'%f',1);scale_F = fscanf(f,'%f',1);N_mat = fscanf(f,'%d',1);for i = 1:N_matn_mat(i) = fscanf(f,'%d',1);ro(i) = fscanf(f,'%f',1);E_el(i) = fscanf(f,'%f',1);mu(i) = fscanf(f,'%f',1);endN_el_type = fscanf(f,'%d',1);el_type = fscanf(f,'%d',N_el_type);rot_vector = fscanf(f,'%d',3);178for i = 1:N_nodesnode_num(i) = fscanf(f,'%d',1);cord(i,:) = fscanf(f,'%e',3);endfor i = 1:N_eln_type_el(i,1) = fscanf(f,'%d',1);n_mat_el(i,1) = fscanf(f,'%d',1);some_num(i,:) = fscanf(f,'%d',8);el_number(i) = fscanf(f,'%d',1);ind_el(i,:) = fscanf(f,'%d',10);endn_fixed_nodes = fscanf(f, '%d', 1);for i = 1:n_fixed_nodesfixations(i,1) = fscanf(f, '%d', 1);endfclose(f);%%%% Executing PSE %%%%%%%% PSE settingsit_number = it_i;format_num = 3;N_S_el_1 = 3000;N_S_el_older = 5000;N_B_S_el_older = 2000;n_freq_desired = n_freq_all;n_precision = 0.00001;n_block = 15;%%% Writing settings file %%%%f_set = fopen([dir_PSE 'settings_eig.txt'], 'w');fprintf(f_set, dir_model, '%s');fprintf(f_set, '\n');fprintf(f_set, ['\\ans_model_' num2str(it_i-1) '.txt'], '%s');fprintf(f_set, '\n');fprintf(f_set, '%d\n', it_number);fprintf(f_set, '%d\n', format_num);fprintf(f_set, '%d\n', N_S_el_1);fprintf(f_set, '%d\n', N_S_el_older);fprintf(f_set, '%d\n', N_B_S_el_older);fprintf(f_set, '%d\n', n_freq_desired);fprintf(f_set, '%f\n', n_precision);fprintf(f_set, '%d\n', n_block);fclose(f_set);%%%% Finding eigenvectors and eigenfrequencies %%%%disp('Calculation started')system('PSE_Client.exe');% clc;disp('Calculation is over')end%%%% Export FE model %%%%function PSE_model_export(f)179global n_DOF N_nodes N_el scale_u scale_F n_mat ro...E_el mu N_el_type el_type rot_vector node_num...cord n_type_el n_mat_el some_num el_number ind_el...n_fixed_nodes fixations N_clustersfprintf(f,'%d\n',n_DOF);fprintf(f,'%d\n',N_nodes);fprintf(f,'%d\n',N_el);fprintf(f,'%f\n',scale_u);fprintf(f,'%f\n',scale_F);fprintf(f,'%d\n',N_clusters);for i = 1:N_clustersfprintf(f,'%d\n',n_mat(i));fprintf(f,'%f\n',ro(i));fprintf(f,'%f\n',E_el(i));fprintf(f,'%f\n',mu(i));endfprintf(f,'%d\n',N_el_type);fprintf(f,'%d\n',el_type);fprintf(f,'%d\n',rot_vector);for i = 1:N_nodesfprintf(f,'%d ',node_num(i));fprintf(f,'%e ',cord(i,:));fprintf(f, '\n');endfor i = 1:N_elfprintf(f,'%d ', n_type_el(i,1));fprintf(f,'%d ', n_mat_el(i,1));fprintf(f,'%d ', some_num(i,:));fprintf(f,'%d ', el_number(i));fprintf(f,'%d ', ind_el(i,:));fprintf(f, '\n');endfprintf(f, '%d\n', n_fixed_nodes);fprintf(f, '%d\n', fixations);fclose(f);end%%%% Results read %%%%function PSE_results_read(f,it_i)global N_el n_freq_all S_EN norm_S_EN freq_i S_EN_0N_el = fread(f, 1, 'int');n_freq = fread(f, 1, 'int');for i = 1:n_freq_allfreq_i(i,1) = fread(f, 1, 'double');modal_mass(i,1) = fread(f, 1, 'double');if (it_i == 1)S_EN_0(:,i) = fread(f, N_el, 'double');else180S_EN_junk = fread(f, N_el, 'double');endendfreq_i = freq_i./(2*pi);fclose(f);end181Приложение С.

Текст программы в среде Matlab для идентификациимодальныхпараметров.Программадляанализаточностистохастической идентификации подпространств SSI EMA_TEST_SSIfunction EMA_TEST_SSIclear all;close all;clc;dt = 0.001;t = 0:dt:5;Fs = 1/dt;ksi1ksi2ksi3ksi4ksi5=====1e-2;2e-2;3e-2;4e-2;5e-2;p1p2p3p4p5=====11;25;33;43;58;A1A2A3A4A5=====1;2;3;4;5;ksi_theory = [ksi1; ksi2; ksi3; ksi4; ksi5];p_theory = [p1; p2; p3; p4; p5];A_theory = [A1; A2; A3; A4; A5];y1y2y3y4y5=====A1*exp(-ksi1*2*pi*p1*t).*sin(2*pi*p1*t);A2*exp(-ksi2*2*pi*p2*t).*sin(2*pi*p2*t);A3*exp(-ksi3*2*pi*p3*t).*sin(2*pi*p3*t);A4*exp(-ksi4*2*pi*p4*t).*sin(2*pi*p4*t);A5*exp(-ksi5*2*pi*p5*t).*sin(2*pi*p5*t);sigma = (0.01*A5)^2;Err = sigma*randn(length(y1),1)';y = y1 + y2 + y3 + y4 + y5 + Err;figure(1)hold on;plot(t,y)xlabel('time')xlim([0 5])L = length(t);NFFT = 2^nextpow2(L);Y = fft(y,NFFT)/L;f = Fs/2*linspace(0,1,NFFT/2+1);A1 = 2*abs(Y(1:NFFT/2+1));алгоритма182figure(2)hold on;plot(f,A1)xlabel('Frequency (Hz)')ylabel('|Y(f)|')xlim([0 80])figure(4)grid on; hold on;plot(p_theory, ksi_theory , 'bo--', 'LineWidth', 1 );%%%%%% SSI %%%%%%%%%shift = 1200;n = 5;[A,A2,C] = SSI_IM_PC_1(y, shift, n);lambda = eig(A);lambda2 = eig(A2);clear A A2mu = log(lambda)./dt;mu2 = log(lambda2)./dt;if (real(mu)>0 & real(mu2)>0)mu = [];mu2 = [];endomega = abs(mu)./(2*pi);ksi = -(real(mu./(2*pi))./omega);omega2 = sqrt(real(mu2).^2+imag(mu2).^2)./(2*pi);ksi2 = -(real(mu2)./omega2);%%%% Sort omega with the accordance to damping %%%%for j = 1:length(omega)for i = 1:length(omega) - jif (omega(i) > omega(i+1))c = omega(i);omega(i) = omega(i+1);omega(i+1) = c;k = ksi(i);ksi(i) = ksi(i+1);ksi(i+1) = k;k = 0;c = 0;endendendksi = ksi(1:2:end-1);p = omega(1:2:end-1);figure(4)hold on; grid on;plot(p,ksi,'ro-')par = [p ksi]183err_p = (abs(p - p_theory)./p_theory).*100;err_ksi = (abs(ksi - ksi_theory)./ksi_theory).*100;err = [err_p err_ksi]end%%%%% SSI procedure %%%%%function[A,A2,C] = SSI_IM_PC_1(Y,i, n)% Y - отклик системы, i сдвиг данных, n - порядок системы %%tic% l - количество датчиков %%[l,N] = size(Y);if(l>N) Y = Y'; [l,N] = size(Y); end% количество столбцов матрицы Ганкеля (j -> inf) %%j = N - 2*i + 1;%%% формирование блочных матриц прошлого и будущего %%Hp = zeros(l*(i-1),j);Hf = Hp;Yp = Hp;%%% Нормировка вэргодичного сигнала %Yn = Y/sqrt(j);clear Yсоответствииопределенияковариационной%%% Заполнение блочных матриц %%for k = 1:iYp((i-k)*l+1:(i-k+1)*l,:) = Yn(:,k:k+j-1);Hf((k-1)*l+1:k*l,:) = Yn(:,i+k:i+k+j-1);Hp((k-1)*l+1:k*l,:) = Yp((i-k)*l+1:(i-k+1)*l,:);endclear Yp;%%% Блочная матрица Теплица %%Ci = Hf*Hp';clear Hp HfCH = Ci*Ci';%%% разложение на собственные значения %%[P1,D1,P1] = svd(CH,0);D = diag(D1);clear D1матрицыдля184%%% Определение порядка системы %%figure(1); hold off;semilogy(D,'b.');title('Singular Values');xlabel('Order');n = 0;while (n < 1) | (n > l*i-1)n = input(' System order ? ');endn = 2*n;%%% Удаление лишних компонент матрицы P %%P1 = P1(:,1:n);%%% Определение матрицы наблюдаемости %%Tmi = P1*diag(sqrt(D(1:n)));%%% Определение динамической матрицы %%Tmi1 = Tmi(1:l*(i-1),:);Tmi2 = Tmi(l+1:end,:);C = Tmi(1:l,:);A = pinv(Tmi1)*Tmi2;A2 = inv(Tmi2'*Tmi1)*Tmi2'*Tmi2;clear Tmi Tmi1 Tmi2tocendОТЗЫВнаучного руководителя соискателя Николаева Сергея Михайловичадоктора технических наук, профессора ВОРОНОВА С.А.Николаев Сергей Михайлович, 1989 г.р., в 2012 г.

Характеристики

Список файлов диссертации