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 = 100e3; % [Pa] ambient pressure, used for load force computation fld.p_cav = 30e3; % [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 = 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_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 = 0; % [-] 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 = NaN; % [Pa] constant in Dowson Higginson relation fld.dow_C2 = NaN; % [-] coefficient in Dowson Higginson relation end % Dynamic viscosity: fld.mu_0 = 3.0e-2; % [Pas] lubricant dynamic viscosity at ambient pressure fld.mu_flag = 0; % [-] 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 = NaN; % [1/Pa] pressure viscosity coefficient, for mineral oils between 1e-8 and 2e-8 [1/Pa] fld.roe_p0 = NaN; % [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 = autorun.sld.flag_elastic(autorun.sim.it); % [-] flag whether solid body is elastic on macroscopic scale (true) or rigid (false) if sld.flag_elastic sld.nu_low = 0.3; % [-] Poisson ratio of lower body sld.E_low = 50e8; % [Pa] Young's modulus of lower body sld.nu_up = 0.3; % [-] Poisson ratio of upper body sld.E_up = 50e8; % [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 end %% Operating condition: opc.flag_imposed = 1; % [-] 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 = 0e-0; % [m] imposed rigid body displacement between upper and lower surface elseif opc.flag_imposed == 2 opc.F_N = NaN; % [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_up_ini = 5; % [m/s] velocity of upper body in x1-direction, the disk opc.u_low_ini = 0; % [m/s] velocity of lower body in x1-direction, the ball [opc.u_up, opc.u_low] = ndgrid(opc.u_up_ini,opc.u_low_ini); % constructing matrices that correspond to each combination of the given initial value opc.u_r = opc.u_up - opc.u_low; % [m/s] relative velocities for each operating condition opc.u_m = (opc.u_up + opc.u_low)/2; % [m/s] mean velocities for each operating condition, must be larger than 0 opc.srr = (opc.u_low - opc.u_up)/(opc.u_low + opc.u_up); % [-] slide-to-roll ratio // -1 for u_low=0 // 0 for u_low=u_up opc.N = numel(opc.u_m); % [-] number of of operating conditions %% Geometry settings: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Study C: Multiple dimples - Woloszynski et al., 2015: % Case a: % Input: % Gap height: geo.h_min = 15e-6; % [m] minimum gap height geo.h_p = 12e-6; % [m] dimple depth % x1-direction: geo.L_x1 = 80e-3; % [m] domain lenght geo.K_x1 = autorun.geo.K_x1(autorun.sim.it); % [-] number of dimples N_unit_x1 = 30; % [-] number of discretization cell for one dimple including surroundings ind_x1_unit_W = 4; % [-] cell index just before dimple ind_x1_unit_E = 27; % [-] cell index just after dimple ind_x1_base_W = 5; % [-] cell index just before dimple base ind_x1_base_E = 26; % [-] cell index just after dimple base % x2-direction: geo.L_x2 = geo.L_x1; geo.K_x2 = geo.K_x1; N_unit_x2 = N_unit_x1; ind_x2_unit_S = ind_x1_unit_W; ind_x2_unit_N = ind_x1_unit_E; ind_x2_base_S = ind_x1_base_W; ind_x2_base_N = ind_x1_base_E; % Computation: N_x1 = geo.K_x1*N_unit_x1; ind_x1 = 1:N_x1; cond_x1_out = mod(ind_x1-1,N_unit_x1)+1 <= ind_x1_unit_W | mod(ind_x1-1,N_unit_x1)+1 >= ind_x1_unit_E; cond_x1_base = mod(ind_x1-1,N_unit_x1)+1 >= ind_x1_base_W+1 & mod(ind_x1-1,N_unit_x1)+1 <= ind_x1_base_E-1; N_x2 = geo.K_x2*N_unit_x2; ind_x2 = 1:N_x2; cond_x2_out = mod(ind_x2-1,N_unit_x2)+1 <= ind_x2_unit_S | mod(ind_x2-1,N_unit_x2)+1 >= ind_x2_unit_N; cond_x2_base = mod(ind_x2-1,N_unit_x2)+1 >= ind_x2_base_S+1 & mod(ind_x2-1,N_unit_x2)+1 <= ind_x2_base_N-1; [cond_x1_out_matr,cond_x2_out_matr] = ndgrid(cond_x1_out,cond_x2_out); [cond_x1_ins_matr,cond_x2_ins_matr] = ndgrid(cond_x1_base,cond_x2_base); clear cond_x1_out;clear cond_x2_out; clear cond_x1_ins; clear cond_x2_ins; cond_out = cond_x1_out_matr==1 | cond_x2_out_matr==1; cond_base = cond_x1_ins_matr==1 & cond_x2_ins_matr==1; clear cond_x1_out_matr;clear cond_x2_out_matr; clear cond_x1_ins_matr; clear cond_x2_ins_matr; cond_par = cond_out | cond_base; cond_incl = ~cond_par; geo.h_g_ma = zeros(N_x1,N_x2); geo.h_g_ma(cond_out==1) = geo.h_min; geo.h_g_ma(cond_base==1) = geo.h_min + geo.h_p; geo.h_g_ma(cond_incl) = geo.h_min + geo.h_p/2; % Add extra points for boundary conditions: h_bound = geo.h_min*ones(N_x1+2,N_x2+2); h_bound(2:N_x1+1,2:N_x2+1) = geo.h_g_ma; geo.h_g_ma = h_bound; clear h_bound; geo.Nx1 = N_x1+2; clear N_x1; geo.dx1 = geo.L_x1/(geo.Nx1 - 1); geo.x1 = linspace(0,geo.L_x1,geo.Nx1); geo.Nx2 = N_x2+2; clear N_x2; geo.dx2 = geo.L_x2/(geo.Nx2 - 1); geo.x2 = linspace(0,geo.L_x2,geo.Nx2); % Save profile information: geo.x3 = -geo.h_g_ma; % [m] macroscopic profile geo.N = geo.Nx1*geo.Nx2; %% Unsteady geometry variation: % Transient algorithm settings: alg.flag_unsteady = false; % [-] boolean, whether simulation is unsteady (true) or steady (false) if alg.flag_unsteady %Define micro-texture parameters: 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 = NaN; % [-] dimensionless time step size, same size as dimensionless space resolution same as Mourier (2006) opc.delta_t = NaN; % [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 = 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 h_t_ma(:,:,it_time,it_OC) = NaN; % [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 = autorun.FBNS.alpha_p_nd(autorun.sim.it); % [-] relaxation factor for update of hydrodynmic pressure alg.FBNS.alpha_thet = autorun.FBNS.alpha_thet(autorun.sim.it); % [-] 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 = NaN; % [-] tolerance of during computation of inital pressure distribution alg.FBNS.p_hd_ini.it_max = NaN; % [-] maximum number of allowed iteration steps before computation of inital pressure distribution is stopped end % Reference quantities: alg.FBNS.ref.x1 = geo.L_x1; % [m] alg.FBNS.ref.x2 = geo.L_x2; % [m] alg.FBNS.ref.h = geo.h_min; % [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 = 1e6; % [Pa] if alg.flag_unsteady alg.FBNS.ref.t = NaN; % [s] end alg.load.toll = alg.FBNS.toll; % [-] tolerance of load balance eqaution if opc.flag_imposed == 2 alg.load.ref.F_N = NaN; % [N] alg.load.h_d_ma_ini = NaN; % [m] initial guess for rigid body displacement alg.load.K_P = NaN; % [-] Propotional coefficient of load balance PID-controller alg.load.K_I = NaN; % [-] Integration coefficient of load balance PID-controller alg.load.K_D = NaN; % [-] 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 = NaN; % [-] 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 fprintf('\n--------------------------------'); fprintf('\nSetup finished!'); fprintf('\nSimulation ID: ''%s''',alg.Sim_ID); fprintf('\n--------------------------------\n');