Add files via upload

Analytical_solutions_DSWIS.m

@@ -0,0 +1,146 @@
+clearvars;
+close all;
+
+%% Analytical solutions for DSWIS (fully symbolic)
+% by Ekaterina Bolotskaya
+
+% 01/15/2022
+
+% This script solves the equation of motion for 1D spring-slider:
+% a = -1/M*(K*(D-V_0*t)+mu*sigma_n)
+% where a - acceleration (second derivative of slip (D))
+% M - mass
+% K - spring stiffness
+% D - slip
+% V_0 - load point velocity
+% t - time
+% mu - friction coefficient
+% sigma_n - normal stress
+% mu is determined by the failure law. In this case double slip-weakening 
+% with initial strengthening (DSWIS) failure law, which consists of one
+% linear strengthening and two linear weakening segments
+
+%% Analytical solution for strengthening (1st segment)
+syms us(t) 
+Dus = diff(us);
+
+syms M K V_0 d_tau_s D_s t real positive
+% d_tau_s = Sn*(mu_s-mu_i) - strengthening stress increase
+% D_s - slip strengthening distance
+
+% Equation
+ode = diff(us,t,2)*M == -K*(us - V_0*t) - d_tau_s*(us)/D_s;
+% Initial conditions
+cond1 = us(0) == 0;
+cond2 = Dus(0) == 0;
+
+% Solve
+uSol(t) = dsolve(ode,[cond1 cond2]);
+
+% Display slip
+uSol = simplify(uSol, 'Steps', 20);
+fprintf('Strengthening:\n u(t) = \n\n');
+pretty(uSol)
+% Display slip rate
+vSol = diff(uSol);
+fprintf('v(t) = \n\n');
+pretty(vSol)
+% Display acceleration
+aSol = diff(vSol);
+fprintf('a(t) = \n\n');
+pretty(aSol)
+
+%% Analytical solution for weakening 1 (2nd segment)
+syms uw1(t) 
+Duw1 = diff(uw1);
+
+syms d_tau_w1 D_w1 real positive
+syms t_s v_s real positive
+% d_tau_w1 = Sn*(mu_s-mu_t) - first weakening stress drop
+% D_w1 - first slip-weakening distance
+% v_s - slip rate at the end of the strengthening segment
+% t_s - time at the end of the strengthening segment
+
+% Equation
+ode = diff(uw1,t,2)*M == -K*(uw1 - V_0*t - V_0*t_s) + d_tau_w1*(uw1)/D_w1 - d_tau_s - K*D_s;
+% Initial conditions
+cond1 = uw1(0) == 0;
+cond2 = Duw1(0) == v_s;
+
+% Solve
+uSol(t) = dsolve(ode,[cond1 cond2]);
+
+% Display slip
+uSol = simplify(uSol, 'Steps', 20);
+fprintf('Strengthening:\n u(t) = \n\n');
+pretty(uSol)
+% Display slip rate
+vSol = diff(uSol);
+fprintf('v(t) = \n\n');
+pretty(vSol)
+% Display acceleration
+aSol = diff(vSol);
+fprintf('a(t) = \n\n');
+pretty(aSol)
+
+%% Analytical solution for weakening 2 (3rd segment)
+syms uw2(t) 
+Duw2 = diff(uw2);
+
+syms d_tau_w2 D_w2 real positive
+syms t_w1 v_w1 real positive
+% d_tau_w2 = Sn*(mu_t-mu_d) - second weakening stress drop
+% D_w2 - second slip-weakening distance
+% v_w1 - slip rate at the end of the first weakening segment
+% t_w1 - time at the end of the first weakening segment
+
+% Equation
+ode = diff(uw2,t,2)*M == -K*(uw2 - V_0*t - V_0*(t_s+t_w1)) + d_tau_w2*(uw2)/D_w2 - d_tau_s+d_tau_w1 - K*(D_s+D_w1);
+% Initial conditions
+cond1 = uw2(0) == 0;
+cond2 = Duw2(0) == v_w1;
+
+% Solve
+uSol(t) = dsolve(ode,[cond1 cond2]);
+
+% Display slip
+uSol = simplify(uSol, 'Steps', 20);
+fprintf('Strengthening:\n u(t) = \n\n');
+pretty(uSol)
+% Display slip rate
+vSol = diff(uSol);
+fprintf('v(t) = \n\n');
+pretty(vSol)
+% Display acceleration
+aSol = diff(vSol);
+fprintf('a(t) = \n\n');
+pretty(aSol)
+
+%% Analytical solution for flat part
+syms uf(t)
+Duf = diff(uf);
+syms t_w2 v_w2 real positive
+% v_w2 - slip rate at the end of the second weakening segment
+% t_w2 - time at the end of the second weakening segment
+
+% Equation
+ode = diff(uf,t,2)*M == -K*(uf - V_0*t - V_0*(t_s+t_w1+t_w2)) - d_tau_s+d_tau_w1+d_tau_w2 - K*(D_s+D_w1+D_w2);
+% Initial conditions
+cond1 = uf(0) == 0;
+cond2 = Duf(0) == v_w2;
+
+% Solve
+uSol(t) = dsolve(ode,[cond1 cond2]);
+
+% Display slip
+uSol = simplify(uSol, 'Steps', 20);
+fprintf('Strengthening:\n u(t) = \n\n');
+pretty(uSol)
+% Display slip rate
+vSol = diff(uSol);
+fprintf('v(t) = \n\n');
+pretty(vSol)
+% Display acceleration
+aSol = diff(vSol);
+fprintf('a(t) = \n\n');
+pretty(aSol)

Analytical_solutions_generic_eqn_linear_friction.m

@@ -0,0 +1,82 @@
+clearvars;
+close all;
+
+%% Analytical solutions for the generic equation (fully symbolic)
+% by Ekaterina Bolotskaya
+
+% 12/30/2021
+
+% This script solves the nondim equation of motion for 1D spring-slider:
+% y" + (1-K_f/K)*y = V_0*t + d_tau_i
+% where y - nondim slip
+% y"      - nondim acceleration
+% t       - nondim time
+% K_f     - failure law segment slope
+% K       - elastic loading slope (spring stiffness)
+% V_0     - nondim load point velocity
+% d_tau_i - nondim stress difference at the beginning of the segment
+% Substitute Ke = (K-K_f)/K - effective loading slope
+
+% There are three "stability regimes" 
+% and three solution types depending on the sign of Ke
+
+%% Equation
+syms y(t) 
+Dy    = diff(y);
+
+syms Ke real
+syms V_0 d_tau_i V_i real 
+% V_i - nondim slip rate at t = 0
+
+ode   = diff(y,t,2) + Ke*y == V_0*t + d_tau_i;
+cond1 = y(0) == 0;
+cond2 = Dy(0) == V_i;
+
+%% Ke > 0
+assumeAlso(Ke > 0)
+% Solve
+ySol(t) = dsolve(ode,[cond1 cond2]);
+
+% Display
+fprintf('K_seg < K \n\n y(t) = \n\n');
+pretty(simplify(ySol))
+
+fprintf('y''(t) = \n\n');
+pretty(simplify(diff(ySol)))
+
+fprintf('y"(t) = \n\n');
+pretty(simplify(diff(diff(ySol))))
+
+%% Ke < 0
+syms Ke real 
+assumeAlso(Ke < 0)
+% Solve
+ySol(t) = dsolve(ode,[cond1 cond2]);
+
+% Display
+fprintf('K_seg > K \n\n y(t) = \n\n');
+pretty(simplify(ySol))
+
+fprintf('y''(t) = \n\n');
+pretty(simplify(diff(ySol)))
+
+fprintf('y"(t) = \n\n');
+pretty(simplify(diff(diff(ySol))))
+
+%% Ke = 0
+% New equation
+ode     = diff(y,t,2) == V_0*t + d_tau_i;
+
+% Solve
+ySol(t) = dsolve(ode,[cond1 cond2]);
+
+% Display
+fprintf('K_seg = K \n\n y(t) = \n\n');
+pretty(simplify(ySol))
+
+fprintf('y''(t) = \n\n');
+pretty(simplify(diff(ySol)))
+
+fprintf('y"(t) = \n\n');
+pretty(simplify(diff(diff(ySol))))
+