classdef Werkstoffmodell < handle properties % model parameters E_A; % elastic modulus of Austenite [N/mm^2]=[GPa] E_M; % elastic Modulus of Martensite [N/mm^2]=[GPa] E_T; % tangent modulus (cf. DIER08) [N/mm^2]=[GPa] E_n_K; % E_n zur Bestimmung von K -------------------------NEU------------------------- eps_Tr_gr_r; % maximum transformation strain [-] eps_Tr_gr_l; % maximum transformation strain [-] sig_fh; % yield stress (forward transformation) [N/mm^2]=[GPa] sig_fr; % yield stress (reverse transformation) [N/mm^2]=[GPa] sig_0; % prestress [N/mm^2] sig_m; % stress at the end of transformation H; % Clausius-Clapeyron coefficient [MPa/K] K; % plastic modulus [N/mm^2] C; % coefficient for temperature rise due to heat of transformation [K] delta; % heat transfer coefficient (Temperatur-Abklingkoeffizient) [1/s] T_0; % temperature at which the yield stresses were determined experimentally [° C] T_R; % room temperature [° C] % values for algorithm: for previous state t_n_1; % previous time [s] eps_n_1_r_stern; % previous overall strain [-] right spring eps_n_1_l_stern; % previous overall strain [-] left spring T_n_1_r; % previous temperature [° C] T_n_1_l; % previous temperature [° C] eps_Tr_n_1_r; % previous transformation strain [-] right eps_Tr_n_1_l; % previous transformation strain [-] left E_n_1_l; % previous elastic modulus [N/mm^2] E_n_1_r; % previous elastic modulus [N/mm^2] sig_n_1_l; % previous stress [N/mm^2] left spring sig_n_1_r; % previous stress [N/mm^2] right spring % values for algorithm: for current state t_n; % current time [s] T_n_r; % current temperature [° C] T_n_l; % current temperature [° C] eps_n_l_stern; % current overall strain [-] left spring eps_n_r_stern; % current overall strain [-] right spring eps_Tr_n_r; % current transformation strain [-] right eps_Tr_n_l; % current transformation strain [-] left E_n_r; % current elastic modulus [N/mm^2] E_n_l; % current elastic modulus [N/mm^2] q_n_r; % current back stress [N/mm^2] q_n_l; % current back stress [N/mm^2] sig_t_n_l; % current test stress [N/mm^2] left spring sig_t_n_r; % current test stress [N/mm^2] right spring sig_n_l_stern; % current stress [N/mm^2] left spring sig_n_r_stern; % current stress [N/mm^2] right spring f_th_n_l; % current transformation function (forward transformation) [N/mm^2] f_th_n_r; % current transformation function (forward transformation) [N/mm^2] f_tr_n_l; % current transformation function (reverse transformation) [N/mm^2] f_tr_n_r; % current transformation function (reverse transformation) [N/mm^2] % values for algorithm: delta values delta_t_n; % difference between the current and previous time [s] delta_T_Tr_n_r; % temperature change due to heat of transformation[°C] delta_T_Tr_n_l; % temperature change due to heat of transformation[°C] delta_eps_n_r; % difference between the current and previous strain [-] delta_eps_n_l; % difference between the current and previous strain [-] delta_eps_Tr_n_r; % change of transformation strain [-] delta_eps_Tr_n_l; % change of transformation strain [-] % not needed for model A_s; % austenite start temperature [K] A_f; % austenite finish temperature [K] M_s; % martensite start temperature [K] M_f; % martensite finish temperature [K] end methods % constructor function obj = Werkstoffmodell(argE_A,argE_M,argE_T,argE_n_K,argKprob,arg_eps_Tr_gr_r,arg_eps_Tr_gr_l,arg_sig_fh,arg_sig_fr,argH,argC,arg_delta,argT_0,argT_R,argA_s,argA_f,argM_s,argM_f) obj.E_A = argE_A; % [N/mm^2]=[MPa] obj.E_M = argE_M; % [N/mm^2]=[MPa] obj.E_T = argE_T; % durch Ablesen; attention when calculating eps_n_1 (maybe: divison by 0) -------------------------NEU------------------------- obj.E_n_K = argE_n_K; % E_n zur Bestimmung von K -------------------------NEU------------------------- obj.K = argKprob; % -------------------------NEU------------------------- obj.eps_Tr_gr_r = arg_eps_Tr_gr_r; % maximum transformation strain [-] obj.eps_Tr_gr_l = arg_eps_Tr_gr_l; % maximum transformation strain [-] obj.sig_fh = arg_sig_fh; % [N/mm^2]=[MPa] obj.sig_fr = arg_sig_fr; % [N/mm^2]=[MPa] obj.H = argH; % Spannungsanstieg durch Temperaturerhöhung (Clausius-Clapeyron-Faktor) obj.C = argC; % "frei waehlbarer Parameter" Achtung: Kelvin oder Grad Celsius? obj.delta = arg_delta; obj.T_0 = argT_0; obj.T_R = argT_R; obj.A_s = argA_s; obj.A_f = argA_f; obj.M_s = argM_s; obj.M_f = argM_f; % initialize the beginning state (standard values) % method fnInitializePreviousState is not needed % eps_n_1 and eps_Tr_n_1 are set with fnSetPrestress % T_n_1 is equal to room temperature and/or can be set with % fnSetRoomTemperature obj.f_th_n_l = 0; obj.f_th_n_r = 0; obj.f_tr_n_l = 0; obj.f_tr_n_r = 0; obj.t_n_1 = 0; obj.eps_n_1_r_stern = 0; obj.eps_n_1_l_stern = 0; obj.T_n_1_r = obj.T_R; obj.T_n_1_l = obj.T_R; obj.eps_Tr_n_1_l = 0; obj.eps_Tr_n_1_r = 0; obj.sig_n_1_l = 0; obj.sig_n_1_r = 0; obj.sig_n_r_stern = 0; obj.sig_n_l_stern = 0; end %% prestress % sets the prestress and calculates the belonging strain and % transformation strain % variables that are influenced: % eps_Tr_n_1 % E_n_1 % eps_n_1 function fnSetPrestress(obj, argPrestress_stern) obj.sig_0 = argPrestress_stern; % prestress at the state n-1 if obj.K ~= 0 obj.sig_m = obj.sig_fh + obj.K * obj.eps_Tr_gr_r; % stress at the end of transformation if obj.sig_0 < obj.sig_fh obj.eps_Tr_n_1_r = 0; obj.eps_Tr_n_1_l = 0; obj.E_n_1_r = obj.E_A; obj.E_n_1_l = obj.E_A; obj.eps_n_1_r_stern = obj.sig_0 / obj.E_n_1_r; obj.eps_n_1_l_stern = obj.sig_0 / obj.E_n_1_l; elseif (obj.sig_0 >= obj.sig_fh) && (obj.sig_0 < obj.sig_m) obj.eps_Tr_n_1_r = (obj.sig_0 - obj.sig_fh) / obj.K; obj.eps_Tr_n_1_l = (obj.sig_0 - obj.sig_fh) / obj.K; obj.E_n_1_r = obj.E_A - obj.eps_Tr_n_1_r / obj.eps_Tr_gr_r * (obj.E_A - obj.E_M); obj.E_n_1_l = obj.E_A - obj.eps_Tr_n_1_l / obj.eps_Tr_gr_l * (obj.E_A - obj.E_M); obj.eps_n_1_r_stern = obj.sig_0 / obj.E_n_1_r + (obj.sig_0 - obj.sig_fh) / obj.K; % --------------???------------------ obj.eps_n_1_l_stern = obj.sig_0 / obj.E_n_1_l + (obj.sig_0 - obj.sig_fh) / obj.K; % --------------???------------------ %obj.eps_n_1_r_stern = obj.sig_fh / obj.E_A + (obj.sig_0 - obj.sig_fh) / obj.K; %obj.eps_n_1_l_stern = obj.sig_fh / obj.E_A + (obj.sig_0 - obj.sig_fh) / obj.K; elseif (obj.sig_0 >= obj.sig_m) obj.eps_Tr_n_1_r = obj.eps_Tr_gr_r; obj.eps_Tr_n_1_l = obj.eps_Tr_gr_l; obj.E_n_1_r = obj.E_M; obj.E_n_1_l = obj.E_M; obj.eps_n_1_r_stern = obj.eps_Tr_gr_r + obj.sig_0 / obj.E_n_1_r; obj.eps_n_1_l_stern = obj.eps_Tr_gr_l + obj.sig_0 / obj.E_n_1_l; end else % K = 0 % Annahme: Ende des Fließens, Fall unrealistisch, perfekte Plastizitaet if obj.sig_0 < obj.sig_fh obj.eps_Tr_n_1_r = 0; obj.eps_Tr_n_1_l = 0; obj.E_n_1_r = obj.E_A; obj.E_n_1_l = obj.E_A; %obj.eps_n_1_r_stern = obj.sig_0 / obj.E_A; %obj.eps_n_1_l_stern = obj.sig_0 / obj.E_A; obj.eps_n_1_r_stern = obj.sig_0 / obj.E_n_1_r; obj.eps_n_1_l_stern = obj.sig_0 / obj.E_n_1_l; elseif obj.sig_0 == obj.sig_fh obj.eps_Tr_n_1_r = obj.eps_Tr_gr_r; obj.eps_Tr_n_1_l = obj.eps_Tr_gr_l; obj.E_n_1_r = obj.E_M; obj.E_n_1_l = obj.E_M; obj.eps_n_1_r_stern = obj.sig_fh / obj.E_M + obj.eps_Tr_gr_r; obj.eps_n_1_l_stern = obj.sig_fh / obj.E_M + obj.eps_Tr_gr_l; elseif obj.sig_0 > obj.sig_fh obj.eps_Tr_n_1_r = obj.eps_Tr_gr_r; obj.eps_Tr_n_1_l = obj.eps_Tr_gr_l; obj.E_n_1_r = obj.E_M; obj.E_n_1_l = obj.E_M; %obj.eps_n_1_r_stern = obj.sig_0 / obj.E_n_1 + obj.eps_Tr_gr_r; %obj.eps_n_1_l_stern = obj.sig_0 / obj.E_n_1 + obj.eps_Tr_gr_l; obj.eps_n_1_r_stern = obj.sig_0 / obj.E_n_1_r + obj.eps_Tr_gr_r; obj.eps_n_1_l_stern = obj.sig_0 / obj.E_n_1_l + obj.eps_Tr_gr_l; end end end %% current state % calculates the current elastic modulus and back stress, as well % as the difference between the current and the previous time and % the difference of the current and previous strain function fnComputeCurrentState(obj, arg_eps_n_r_stern, arg_eps_n_l_stern, arg_t_n) obj.eps_n_r_stern = arg_eps_n_r_stern; obj.eps_n_l_stern = arg_eps_n_l_stern; obj.t_n = arg_t_n; obj.delta_t_n = obj.t_n - obj.t_n_1; obj.delta_eps_n_r = obj.eps_n_r_stern - obj.eps_n_1_r_stern; obj.delta_eps_n_l = obj.eps_n_l_stern - obj.eps_n_1_l_stern; obj.q_n_r = obj.H * (obj.T_n_1_r - obj.T_0); obj.q_n_l = obj.H * (obj.T_n_1_l - obj.T_0); obj.E_n_r = obj.E_A - obj.eps_Tr_n_1_r / obj.eps_Tr_gr_r * (obj.E_A - obj.E_M); obj.E_n_l = obj.E_A - obj.eps_Tr_n_1_l / obj.eps_Tr_gr_l * (obj.E_A - obj.E_M); end %% predictor step % calculates the current stress, the foward transformation % function and the reverse transformation function for the test % case (may be corrected later), respectively function fnRunPredictorStep(obj) obj.sig_t_n_r = obj.E_n_r * (obj.eps_n_r_stern - obj.eps_Tr_n_1_r); % -----------------obj.eps_Tr_n_1_r ??------------------ obj.sig_t_n_l = obj.E_n_l * (obj.eps_n_l_stern - obj.eps_Tr_n_1_l); % -----------------obj.eps_Tr_n_1_l ??------------------ %obj.sig_t_n_r = obj.E_n_r * (obj.eps_n_r_stern - obj.eps_n_1_r_stern); %obj.sig_t_n_l = obj.E_n_l * (obj.eps_n_l_stern - obj.eps_n_1_l_stern); obj.f_th_n_r = obj.sig_t_n_r - (obj.sig_fh + obj.q_n_r + obj.K * obj.eps_Tr_n_1_r); obj.f_th_n_l = obj.sig_t_n_l - (obj.sig_fh + obj.q_n_l + obj.K * obj.eps_Tr_n_1_l); obj.f_tr_n_r = obj.sig_t_n_r - (obj.sig_fr + obj.q_n_r + obj.K * obj.eps_Tr_n_1_r); obj.f_tr_n_l = obj.sig_t_n_l - (obj.sig_fr + obj.q_n_l + obj.K * obj.eps_Tr_n_1_l); end %% flow conditions % evaluates the flow condition [SCHM04, S.92] and runs the % appropriate corrector step function fnEvaluateFlowCondition(obj) %Rechte Seite if (obj.delta_eps_n_r > 0) && (obj.f_th_n_r <= 0) && (obj.eps_Tr_n_1_r <= obj.eps_Tr_gr_r) % Fall a) obj.fnRunElasticCorrectorStep_r; elseif (obj.delta_eps_n_r > 0) && (obj.f_th_n_r > 0) && (obj.eps_Tr_n_1_r < obj.eps_Tr_gr_r) % Fall b) obj.fnRunForwardCorrectorStep_r; elseif (obj.delta_eps_n_r > 0) && (obj.f_th_n_r > 0) && (obj.eps_Tr_n_1_r == obj.eps_Tr_gr_r) % Fall c) obj.fnRunElasticCorrectorStep_r; elseif (obj.delta_eps_n_r < 0) && (obj.f_tr_n_r >= 0) && (obj.eps_Tr_n_1_r >= 0) % Fall d) obj.fnRunElasticCorrectorStep_r; elseif (obj.delta_eps_n_r < 0) && (obj.f_tr_n_r < 0) && (obj.eps_Tr_n_1_r > 0) % Fall e) obj.fnRunReverseCorrectorStep_r; elseif (obj.delta_eps_n_r < 0) && (obj.f_tr_n_r < 0) && (obj.eps_Tr_n_1_r == 0) % Fall f) obj.fnRunElasticCorrectorStep_r; % elseif (obj.delta_eps_n == 0) % Fall benötigt, um Zustand zum Zeitpunkt t = 0s ausgeben zu lassen % obj.fnRunElasticCorrectorStep; end %Linke Seite if (obj.delta_eps_n_l > 0) && (obj.f_th_n_l <= 0) && (obj.eps_Tr_n_1_l <= obj.eps_Tr_gr_l) % Fall a) obj.fnRunElasticCorrectorStep_l; elseif (obj.delta_eps_n_l > 0) && (obj.f_th_n_l > 0) && (obj.eps_Tr_n_1_l < obj.eps_Tr_gr_l) % Fall b) obj.fnRunForwardCorrectorStep_l; elseif (obj.delta_eps_n_l > 0) && (obj.f_th_n_l > 0) && (obj.eps_Tr_n_1_l == obj.eps_Tr_gr_l) % Fall c) obj.fnRunElasticCorrectorStep_l; elseif (obj.delta_eps_n_l < 0) && (obj.f_tr_n_l >= 0) && (obj.eps_Tr_n_1_l >= 0) % Fall d) obj.fnRunElasticCorrectorStep_l; elseif (obj.delta_eps_n_l < 0) && (obj.f_tr_n_l < 0) && (obj.eps_Tr_n_1_l > 0) % Fall e) obj.fnRunReverseCorrectorStep_l; elseif (obj.delta_eps_n_l < 0) && (obj.f_tr_n_l < 0) && (obj.eps_Tr_n_1_l == 0) % Fall f) obj.fnRunElasticCorrectorStep_l; % elseif (obj.delta_eps_n == 0) % Fall benötigt, um Zustand zum Zeitpunkt t = 0s ausgeben zu lassen % obj.fnRunElasticCorrectorStep; end end %% elastic corrector step RIGHT % sets the current stress and transformation strain for the elastic % case function fnRunElasticCorrectorStep_r(obj) obj.eps_Tr_n_r = obj.eps_Tr_n_1_r; obj.sig_n_r_stern = obj.sig_t_n_r; % disp('elastic'); end %% elastic corrector step LEFT % sets the current stress and transformation strain for the elastic % case function fnRunElasticCorrectorStep_l(obj) obj.eps_Tr_n_l = obj.eps_Tr_n_1_l; obj.sig_n_l_stern = obj.sig_t_n_l; % disp('elastic'); end %% corrector step for foward transformation RIGHT % sets the current stress and transformation strain in case of a % forward transformation function fnRunForwardCorrectorStep_r(obj) obj.delta_eps_Tr_n_r = obj.f_th_n_r / (obj.E_n_r + obj.K + obj.H * obj.C); if (0 <= obj.eps_Tr_n_1_r + obj.delta_eps_Tr_n_r) && (obj.eps_Tr_n_1_r + obj.delta_eps_Tr_n_r <= obj.eps_Tr_gr_r) obj.eps_Tr_n_r = obj.eps_Tr_n_1_r + obj.delta_eps_Tr_n_r; elseif (obj.eps_Tr_n_1_r + obj.delta_eps_Tr_n_r < 0) obj.eps_Tr_n_r = 0; elseif (obj.eps_Tr_n_1_r + obj.delta_eps_Tr_n_r > obj.eps_Tr_gr_r) obj.eps_Tr_n_r = obj.eps_Tr_gr_r; end obj.sig_n_r_stern = obj.sig_t_n_r - obj.E_n_r * obj.delta_eps_Tr_n_r; % disp('forward'); end %% corrector step for foward transformation LEFT % sets the current stress and transformation strain in case of a % forward transformation function fnRunForwardCorrectorStep_l(obj) obj.delta_eps_Tr_n_l = obj.f_th_n_l / (obj.E_n_l + obj.K + obj.H * obj.C); if (0 <= obj.eps_Tr_n_1_l + obj.delta_eps_Tr_n_l) && (obj.eps_Tr_n_1_l + obj.delta_eps_Tr_n_l <= obj.eps_Tr_gr_l) obj.eps_Tr_n_l = obj.eps_Tr_n_1_l + obj.delta_eps_Tr_n_l; elseif (obj.eps_Tr_n_1_l + obj.delta_eps_Tr_n_l < 0) obj.eps_Tr_n_l = 0; elseif (obj.eps_Tr_n_1_l + obj.delta_eps_Tr_n_l > obj.eps_Tr_gr_l) obj.eps_Tr_n_l = obj.eps_Tr_gr_l; end obj.sig_n_l_stern = obj.sig_t_n_l - obj.E_n_l * obj.delta_eps_Tr_n_l; % disp('forward'); end %% corrector step for reverse transformation RIGHT % sets the current stress and transformation strain in case of a % reverse transformation function fnRunReverseCorrectorStep_r(obj) obj.delta_eps_Tr_n_r = obj.f_tr_n_r / (obj.E_n_r + obj.K + obj.H * obj.C); if (0 <= obj.eps_Tr_n_1_r + obj.delta_eps_Tr_n_r) && (obj.eps_Tr_n_1_r + obj.delta_eps_Tr_n_r <= obj.eps_Tr_gr_r) obj.eps_Tr_n_r = obj.eps_Tr_n_1_r + obj.delta_eps_Tr_n_r; elseif (obj.eps_Tr_n_1_r + obj.delta_eps_Tr_n_r < 0) obj.eps_Tr_n_r = 0; elseif (obj.eps_Tr_n_1_r + obj.delta_eps_Tr_n_r > obj.eps_Tr_gr_r) obj.eps_Tr_n_r = obj.eps_Tr_gr_r; end obj.sig_n_r_stern = obj.sig_t_n_r - obj.E_n_r * obj.delta_eps_Tr_n_r; % disp('reverse'); end %% corrector step for reverse transformation LEFT % sets the current stress and transformation strain in case of a % reverse transformation function fnRunReverseCorrectorStep_l(obj) obj.delta_eps_Tr_n_l = obj.f_tr_n_l / (obj.E_n_l + obj.K + obj.H * obj.C); if (0 <= obj.eps_Tr_n_1_l + obj.delta_eps_Tr_n_l) && (obj.eps_Tr_n_1_l + obj.delta_eps_Tr_n_l <= obj.eps_Tr_gr_l) obj.eps_Tr_n_l = obj.eps_Tr_n_1_l + obj.delta_eps_Tr_n_l; elseif (obj.eps_Tr_n_1_l + obj.delta_eps_Tr_n_l < 0) obj.eps_Tr_n_l = 0; elseif (obj.eps_Tr_n_1_l + obj.delta_eps_Tr_n_l > obj.eps_Tr_gr_l) obj.eps_Tr_n_l = obj.eps_Tr_gr_l; end obj.sig_n_l_stern = obj.sig_t_n_l - obj.E_n_l * obj.delta_eps_Tr_n_l; % disp('obj.eps_n_l_stern ='); disp(obj.eps_n_l_stern); % disp('obj.sig_n_l_stern ='); disp(obj.sig_n_l_stern); % % disp('!!! test !!!'); % disp('!!! test !!!'); % disp('reverse'); end %% update temperature % sets the current temperature and the temperature difference due % to transformation in contrast to the previous step function fnUpdateTemperature(obj) obj.delta_T_Tr_n_r = obj.C * (obj.eps_Tr_n_r - obj.eps_Tr_n_1_r); obj.delta_T_Tr_n_l = obj.C * (obj.eps_Tr_n_l - obj.eps_Tr_n_1_l); obj.T_n_r = (obj.T_n_1_r + obj.delta_T_Tr_n_r - obj.T_R) * exp(-obj.delta * obj.delta_t_n) + obj.T_R; obj.T_n_l = (obj.T_n_1_l + obj.delta_T_Tr_n_l - obj.T_R) * exp(-obj.delta * obj.delta_t_n) + obj.T_R; end %% update state variables % saves the state variables for time, strain, temperature and % transformation strain, respectively function fnUpdateStateVariables(obj) obj.t_n_1 = obj.t_n; obj.eps_n_1_r_stern = obj.eps_n_r_stern; obj.eps_n_1_l_stern = obj.eps_n_l_stern; obj.T_n_1_r = obj.T_n_r; obj.T_n_1_l = obj.T_n_l; obj.eps_Tr_n_1_r = obj.eps_Tr_n_r; obj.eps_Tr_n_1_l = obj.eps_Tr_n_l; obj.sig_n_1_r = obj.sig_n_r_stern; obj.sig_n_1_l = obj.sig_n_l_stern; end end end