close all; clc; clearvars -except autorun; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % EHL setup code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Main developer: % Erik Hansen: erik.hansen@kit.edu % Further developers who contributed to the code: % Altay Kacan, 08.01.2021: Implementation of the unsteady terms in the % Reynolds equation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Information: % This code creates and saves a setup to be then used by the script % EHL_02_mainprocess.m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Simulation Identifier: alg.Sim_ID = autorun.alg.Sim_ID{autorun.sim.it}; %% Fluid settings: fld.h_min = 1e-9; % [m] minimum allowed gap height to prevent it from becoming zero fld.p_amb = 0e0; % [Pa] ambient pressure, used for load force computation fld.p_cav = fld.p_amb; % [Pa] cavitation pressure % Boundary conditions: % Hydrodynamic pressure: fld.flag_p_bound = 0; % [-] flag, which boundary conditions for the hydrodynamic pressure are used: % 0) Dirichlet at West, North, South and East boundary % 1) Dirichlet at West and East boundary, Neumann at North and South boundary % 2) Dirichlet at West boundary, Neumann at North, South and East boundary if fld.flag_p_bound == 0 fld.p_S_side = fld.p_amb; % [Pa] pressure on the South side, must be in interval of [fld.p_cav,inf] fld.p_W_side = fld.p_amb; % [Pa] pressure at the inlet, must be in interval of [fld.p_cav,inf] fld.p_E_side = fld.p_amb; % [Pa] pressure at the outlet, must be in interval of [fld.p_cav,inf] fld.p_N_side = fld.p_amb; % [Pa] pressure on the North side, must be in interval of [fld.p_cav,inf] elseif fld.flag_p_bound == 1 fld.p_S_side = NaN; % [Pa] pressure on the South side, must be in interval of [fld.p_cav,inf] fld.p_W_side = fld.p_amb; % [Pa] pressure at the inlet, must be in interval of [fld.p_cav,inf] fld.p_E_side = fld.p_amb; % [Pa] pressure at the outlet, must be in interval of [fld.p_cav,inf] fld.p_N_side = NaN; % [Pa] pressure on the North side, must be in interval of [fld.p_cav,inf] elseif fld.flag_p_bound == 2 fld.p_S_side = NaN; % [Pa] pressure on the South side, must be in interval of [fld.p_cav,inf] fld.p_W_side = fld.p_amb; % [Pa] pressure at the inlet, must be in interval of [fld.p_cav,inf] fld.p_E_side = NaN; % [Pa] pressure at the outlet, must be in interval of [fld.p_cav,inf] fld.p_N_side = NaN; % [Pa] pressure on the North side, must be in interval of [fld.p_cav,inf] end % Cavity fraction: fld.flag_thet_bound = 2; % [-] flag, which boundary conditions for the hydrodynamic pressure are used: % 0) Dirichlet at West, North, South and East boundary % 1) Dirichlet at West and East boundary, Neumann at North and South boundary % 2) Dirichlet at West boundary, Neumann at North, South and East boundary if fld.flag_thet_bound == 0 fld.thet_S_side = 0.0; % [-] cavity fraction on the South side, must be in interval of [0,1] fld.thet_W_side = 0.0; % [-] cavity fraction at the inlet, must be in interval of [0,1] fld.thet_E_side = 0.0; % [-] cavity fraction at the outlet, must be in interval of [0,1] fld.thet_N_side = 0.0; % [-] cavity fraction on the North side, must be in interval of [0,1] elseif fld.flag_thet_bound == 1 fld.thet_S_side = NaN; % [-] cavity fraction on the South side, must be in interval of [0,1] fld.thet_W_side = 0.0; % [-] cavity fraction at the inlet, must be in interval of [0,1] fld.thet_E_side = 0.0; % [-] cavity fraction at the outlet, must be in interval of [0,1] fld.thet_N_side = NaN; % [-] cavity fraction on the North side, must be in interval of [0,1] elseif fld.flag_thet_bound == 2 fld.thet_S_side = NaN; % [-] cavity fraction on the South side, must be in interval of [0,1] fld.thet_W_side = 0.0; % [-] cavity fraction at the inlet, must be in interval of [0,1] fld.thet_E_side = NaN; % [-] cavity fraction at the outlet, must be in interval of [0,1] fld.thet_N_side = NaN; % [-] cavity fraction on the North side, must be in interval of [0,1] end % Density: fld.rho_0 = 850; % [kg/m^3] lubricant density at ambient pressure, value does not actually matter since dimensionless Reynolds equation is solved fld.rho_flag = 1; % [-] flag whether liquid lubricant density is constant (0) or Dowson-Higginson relation is used for pressure dependency (1) if fld.rho_flag == 1 fld.dow_C1 = 5.9e8; % [Pa] constant in Dowson Higginson relation fld.dow_C2 = 1.34; % [-] coefficient in Dowson Higginson relation end % Dynamic viscosity: fld.mu_0 = 2.5e-1; % [Pas] lubricant dynamic viscosity at ambient pressure fld.mu_flag = 2; % [-] flag whether liquid lubricant dynamic viscosity is constant (0), Barus relation (1) or Roelands relation (2) is used for pressure dependency if fld.mu_flag == 1 fld.barus_alpha = NaN; % [1/Pa] pressure viscosity coefficient, for mineral oils between 1e-8 and 2e-8 [1/Pa] elseif fld.mu_flag == 2 fld.roe_alpha = 22e-9; % [1/Pa] pressure viscosity coefficient, for mineral oils between 1e-8 and 2e-8 [1/Pa] fld.roe_p0 = 1.96e8; % [Pa] constant in Roelands equation fld.roe_z = fld.roe_alpha*fld.roe_p0/(log(fld.mu_0) + 9.67); % [-] pressure viscosity index end % Shear thinning: fld.shear_thinning_model_flag = 0; % [-] flag how shear thinning is considered % 0: no shear thinning % 1: cut-off % 2: Eyring if fld.shear_thinning_model_flag > 0 fld.shear_thinning_limit_flag = 1; % [-] flag how shear thinning limit is considered if fld.shear_thinning_limit_flag == 1 % Constant limit fld.tau_max = NaN; % [Pa] limiting shear stress elseif fld.shear_thinning_limit_flag == 2 % Varying limit fld.tau_max_coeff = NaN; % [-] fraction of hydrodynamic pressure that is used for the limiting shear stress end end %% Solid settings: sld.flag_elastic = true; % [-] flag whether solid body is elastic on macroscopic scale (true) or rigid (false) if sld.flag_elastic sld.nu_low = NaN; % [-] Poisson ratio of lower body sld.E_low = NaN; % [Pa] Young's modulus of lower body sld.nu_up = NaN; % [-] Poisson ratio of upper body sld.E_up = NaN; % [Pa] Young's modulus of upper body % sld.E_dash = 2/((1 - sld.nu_low^2)/sld.E_low + (1 - sld.nu_up^2)/sld.E_up); % [Pa] reduced modulus of elasticity sld.E_dash = 110e9; % [Pa] reduced modulus of elasticity end %% Operating condition: opc.flag_imposed = 2; % [-] flag whether a constant rigid body displacement h_d (1) or a constant normal load force W (2) is imposed if opc.flag_imposed == 1 opc.h_d = NaN; % [m] imposed rigid body displacement between upper and lower surface elseif opc.flag_imposed == 2 opc.F_N = 15; % [N] imposed normal load force in negative x3-direction end % To specify several operating conditions, just set up opc.u_m_ini and/or % opc.u_r_ini as 1D arrays, e.g.: opc.u_m_ini = [9.0e-2 1.0e-2]; opc.u_m_ini = autorun.opc.u_m_ini(autorun.sim.it); % [m/s] mean velocities for each operating condition, must be larger than 0 opc.srr_ini = autorun.opc.srr_ini(autorun.sim.it); % [-] slide-to-roll ratio // -1 for u_low=0 // 0 for u_low=u_up [opc.u_m, opc.srr] = ndgrid(opc.u_m_ini,opc.srr_ini); % constructing matrices that correspond to each combination of the given initial value opc.u_r = -2.*opc.u_m.*opc.srr; % [m/s] relative velocities for each operating condition opc.u_up = opc.u_m + opc.u_r/2; % [m/s] velocity of upper body in x1-direction, the disk opc.u_low = opc.u_m - opc.u_r/2; % [m/s] velocity of lower body in x1-direction, the ball opc.N = numel(opc.u_m); % [-] number of of operating conditions %% Geometry settings: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Study A: Ball on Disc - Mourier et al., 2006: % Define geometry parameters geo.Nx1 = 257; % [-] number of discretization points in the x1-direction geo.Nx2 = 257; % [-] number of discretization points in the x2-direction geo.a = 136.5e-6; % [m] radius of Hertzian contact from the experimental setup of Mourier (2006) x1_Wb = -3*geo.a; % [m] x1-coordinate of cell center at the West boundary x1_Eb = 3*geo.a; % [m] x1-coordinate of cell center at the East boundary x2_Sb = -3*geo.a; % [m] x2-coordinate of cell center at the South boundary x2_Nb = 3*geo.a; % [m] x2-coordinate of cell center at the North boundary geo.Rx1 = 12.5e-3; % [m] ball radius % Construct geometry: geo.dx1 = (x1_Eb - x1_Wb)/(geo.Nx1 - 1); % [m] spacing in x1-direction geo.dx2 = (x2_Nb - x2_Sb)/(geo.Nx2 - 1); % [m] spacing in x2-direction geo.x1 = linspace(x1_Wb,x1_Eb,geo.Nx1); % [m] x1-coordinates with uniform discretization geo.x2 = linspace(x2_Sb,x2_Nb,geo.Nx2); % [m] x2-coordinates with uniform discretization [x1_matr,x2_matr] = ndgrid(geo.x1,geo.x2); geo.h_g_ma = (x1_matr.^2)/(2*geo.Rx1) + (x2_matr.^2)/(2*geo.Rx1); % [m] gap height variation induced by macroscopic profile without the microtexture % Save profile information: geo.x3 = -geo.h_g_ma; % [m] macroscopic profile geo.N = geo.Nx1*geo.Nx2; %% Unsteady geometry variation: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Study A: Ball on Disc - Mourier et al., 2006: % Transient algorithm settings: alg.flag_unsteady = true; % [-] boolean, whether simulation is unsteady (true) or steady (false) if alg.flag_unsteady %Define micro-texture parameters: geo_Ti.r = autorun.geo_Ti.r(autorun.sim.it); % [m] microcavity radius (different for fig 5 & 8 in Mourier (2006)) geo_Ti.d = autorun.geo_Ti.d(autorun.sim.it); % [m] microcavity depth from Mourier (different for simulation & experiment) geo_Ti.x1_0 = geo.x1(1); % [m] initial x1 coordinate of the micro-texture center geo_Ti.x2_0 = 0; % [m] initial x2 coordinate of the micro-texture center geo_Ti.u_tex = opc.u_low; % [m/s] velocity of the micro-texture delta_t_nd = geo.dx1/geo.a; % [-] dimensionless time step size, same size as dimensionless space resolution same as Mourier (2006) opc.delta_t = geo.a./(opc.u_m).*delta_t_nd; % [s] time step size for each operating condition opc.N_t = round((geo.x1(geo.Nx1) - geo.x1(1))./(opc.delta_t.*geo_Ti.u_tex) + 1); % [-] amount of discrete time steps % Initialize the gap height change and the matrix containing the dimple center coordinates: time = zeros(max(opc.N_t(:)),opc.N); % [s] discrete times h_t_ma = zeros(geo.Nx1,geo.Nx2,max(opc.N_t(:)),opc.N); % [m] 4 dimensional array (2 spatial, 1 temporal, 1 extra dimension for operating condition) x1_c_vec = zeros(max(opc.N_t(:)),opc.N); % [m] 2 dimensional array x1 coordinates of the dimple center for each time point and operating condition x2_c_vec = zeros(max(opc.N_t(:)),opc.N); % [m] 2 dimensional array x2 coordinates of the dimple center for each time point and operating condition % Loop for going through all operating conditions: for it_OC = 1:opc.N % Time loop to precalculate imposed transient gap height change for all time points: for it_time = 1:opc.N_t(it_OC) fprintf("Precalculating time point: %i/%i for operating condition: %i/%i\n",it_time,opc.N_t(it_OC),it_OC,opc.N); % Calculate corresponding time: time(it_time,it_OC) = (it_time-1)*opc.delta_t(it_OC); % At the first time point time is = 0, the time step size depends on the operating condition % Calculate gap height change due to micro-texture for this time point: x1_c = geo_Ti.x1_0 + geo_Ti.u_tex(it_OC)*time(it_time,it_OC); % [m] time dependent x1 coordinate of the texture center x1_c_vec(it_time,it_OC) = x1_c; % [m] saved value of the x1 coordinate of the dimple center x2_c = geo_Ti.x2_0; % [m] x2 coordinate of the texture center x2_c_vec(it_time,it_OC) = x2_c; % [m] saved value of the x1 coordinate of the dimple center dist = sqrt((x1_matr-x1_c).^2+(x2_matr-x2_c).^2); % [m] distance of grid points to the center of the micro-texture h_t_ma(:,:,it_time,it_OC) = geo_Ti.d*cos(pi/2*dist/(1.2*geo_Ti.r)).*exp(... -2*(dist/(1.2*geo_Ti.r)).^2); % [m] gap height influence due to the micro-texture, the third dimension is time and the fourth dimension is the corresponding operating condition end end else opc.N_t = ones(size(opc.u_m,1),size(opc.u_m,2));% [-] amount of discrete time steps end %% Algorithm settings: % Genereal relative tolerance: alg.it_tot_max = 1e3; % [-] maximum number of total iterations until simulation is aborted alg.FBNS.toll = 1e-6; % [-] tolerance for FBNS residual checks alg.FBNS.alpha_p_nd = 5e-2; % [-] relaxation factor for update of hydrodynmic pressure alg.FBNS.alpha_thet = 5e-2; % [-] relaxation factor for update of cavity fraction alg.FBNS.flag_enforce_limit_p_nd_red_min = true; % [-] boolean whether to truncate updated p_nd_red_min at 0 (true) or not (false) alg.FBNS.flag_enforce_limit_thet_red_min = true; % [-] boolean whether to truncate updated thet_red_min at 0 (true) or not (false) alg.FBNS.flag_enforce_limit_p_red_max = false; % [-] boolean whether to truncate updated p_nd_red_min at alg.FBNS.limit_p_nd_red_max (true) or not (false) if alg.FBNS.flag_enforce_limit_p_red_max alg.FBNS.limit_p_red_max = NaN; % [Pa] maximum limit for the hydrdynamic pressure end alg.FBNS.flag_enforce_limit_thet_red_max = true; % [-] boolean whether to truncate updated thet_red_min at 1 (true) or not (false) alg.FBNS.flag_res = 10; % [-] flag, which residual definition is used % The FBNS residual defintions can be found and added in the function FBNS_det_residual % When in doubt choose alg.FBNS.flag_res = 10 (strict with respect to, G,F,p_hd and theta) % or alg.FBNS.flag_res = 11 (strict with respect to, G,p_hd) alg.FBNS.flag_print_it = 1; % [-] flag, whether information about each iteration of the FBNS algorithm is printed (1) or not (0) alg.FBNS.flag_Couette_Discrete = autorun.alg.FBNS.flag_Couette_Discrete(autorun.sim.it); % [-] flag, which discretization scheme is used for the Couette term % (maximum order is 2 since midpoint rule is always used for integral flux approximations): % 1) 1.Order Upwind, 2) 2.Order QUICK, 3) 2.Order LUI, 4) 2.Order CUI alg.FBNS.flag_A_h_S = true; % [-] boolean whether (true) or not (false) to use the South Diagonal of A_h: true always advised alg.FBNS.flag_A_h_WW = false; % [-] boolean whether (true) or not (false) to use the WestWest Diagonal of A_h: true advised if LUI is used, otherwise false is advised alg.FBNS.flag_A_h_E = true; % [-] boolean whether (true) or not (false) to use the East Diagonal of A_h: true always advised alg.FBNS.flag_A_h_N = true; % [-] boolean whether (true) or not (false) to use the North Diagonal of A_h: true always advised alg.FBNS.flag_p_hd_ini = 1; % [-] flag, which initial guess is used for the hydrodynamic pressure distribution % 1) fld.p_amb, 2) elastic dry contact solution if alg.FBNS.flag_p_hd_ini == 2 alg.FBNS.p_hd_ini.toll = alg.FBNS.toll; % [-] tolerance of during computation of inital pressure distribution alg.FBNS.p_hd_ini.it_max = 2e2; % [-] maximum number of allowed iteration steps before computation of inital pressure distribution is stopped end % Reference quantities: alg.FBNS.ref.x1 = geo.a; % [m] alg.FBNS.ref.x2 = geo.a; % [m] alg.FBNS.ref.h = ((geo.a)^2)/geo.Rx1; % [m] alg.FBNS.ref.mu = fld.mu_0; % [Pas] alg.FBNS.ref.rho = fld.rho_0; % [kg/m^3] alg.FBNS.ref.u = opc.u_m; % [m/s] alg.FBNS.ref.p = 385e6; % [Pa] if alg.flag_unsteady alg.FBNS.ref.t = geo.a./opc.u_m; % [s] end alg.load.toll = alg.FBNS.toll; % [-] tolerance of load balance eqaution if opc.flag_imposed == 2 alg.load.ref.F_N = opc.F_N; % [N] alg.load.h_d_ma_ini = -1e-6; % [m] initial guess for rigid body displacement alg.load.K_P = 1e-3; % [-] Propotional coefficient of load balance PID-controller alg.load.K_I = 2e-2; % [-] Integration coefficient of load balance PID-controller alg.load.K_D = 1e-3; % [-] Derivative coefficient of load balance PID-controller end % Advanced algorithm settings: alg.flag_profile_viewer = false; % [-] boolean, whether profile viewer is activated (true) or not (false) alg.flag_limit_num_comp_threads = false; % [-] boolean, whether whther the number of computational threads is limited (true) or not (false) if alg.flag_limit_num_comp_threads alg.max_num_comp_threads = 8; % [-] maximum number of computational threads end %% Save setup data: % Define path: output_path = './../data/EHL_02_mainprocess/Input/'; % Create folder/ delete previous contents: warning('off','MATLAB:MKDIR:DirectoryExists') mkdir (output_path) rmdir(output_path, 's') % remove old setups if they exist mkdir (output_path) % Save data: save (fullfile(output_path,'fld.mat'),'fld'); save (fullfile(output_path,'sld.mat'),'sld'); save (fullfile(output_path,'geo.mat'),'geo'); save (fullfile(output_path,'opc.mat'),'opc'); save (fullfile(output_path,'alg.mat'),'alg'); if alg.flag_unsteady save (fullfile(output_path,'geo_Ti.mat'),'geo_Ti'); for it_OC=1:opc.N % Create directories for each operating condition: output_path_OC = fullfile(output_path,sprintf('h_time/OC_%i/',it_OC)); mkdir (output_path_OC) for it_time=1:opc.N_t(it_OC) % Create directories for each time point: output_path_time = fullfile(output_path_OC,sprintf('Time_%i/',it_time)); mkdir (output_path_time) % Create and save the corresponding slices of the 4D array: slc.h_t_ma = h_t_ma(:,:,it_time,it_OC); slc.x1_c = x1_c_vec(it_time,it_OC); slc.x2_c = x2_c_vec(it_time,it_OC); slc.t = time(it_time,it_OC); save(fullfile(output_path_time,'slc.mat'),'slc'); end end end %% Plot settings: flag_plot_setup_data = false; if flag_plot_setup_data % Line plots: KIT_colorlist={[0,150,130]/255,[162 34 35]/255,[70 100 170]/255,[252 229 0]/255,[140 182 60]/256,[223 155 27]/255,[167 130 46]/255,[163 16 124]/255,[35 161 224]/255}; % Surface plots: colmap = parula; % Size: widthlines = 1.2; sizeoffonts = 11; sizeoflegendfonts = 11; % Plots: % Surface plot gap height: fig = figure('Units','centimeters','Position',[10 08 7 7]); if alg.flag_unsteady surf(geo.x1,geo.x2,geo.h_g_ma'+h_t_ma(:,:,1,1)') else surf(geo.x1,geo.x2,geo.h_g_ma') end xlabel('x_1 [m]','fontsize',sizeoffonts); ylabel('x_2 [m]','fontsize',sizeoffonts); zlabel('h_m [m]','fontsize',sizeoffonts); title('Mean gap height','fontsize',sizeoffonts) fig.CurrentAxes.ZDir = 'Reverse'; shading interp colormap(flipud(colmap)) light lighting 'gouraud' lightangle(-90,-30) clear fig; end fprintf('\n--------------------------------'); fprintf('\nSetup finished!'); fprintf('\nSimulation ID: ''%s''',alg.Sim_ID); fprintf('\n--------------------------------\n');