• API
• FAQ
• Tools
• Archive
daily pastebin goal
23%
SHARE
TWEET

# Untitled

a guest Oct 21st, 2017 44 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. % mesh
2.     DL = 0.08 ;
3.     sd_ratio = 3;
4.     SL = DL*sd_ratio ;
5.     Np = 8 ;
6.     N = floor(Np/DL);
7.     dx = L/N;
8.     a = dx ;
9.     x = zeros(1,N+1) ;
10.     q = ones(1, N) ;
11.     G = ones(1, N) ;
12.     lam = 1e-5 ; %
13.     beta = ro/EI ;
14.     dt= sqrt(lam*beta)*dx^2 ;
15.
16.     % Boundaries
17.     yp1(N+1) = 0 ;
18.     y = 0 ;
19.     ym1 = 0 ;
20.
21.     % Implementation
22.     i = 1 : N+1 ;
23.     x = dx * (i-1) ;
24.
25.     tf = 0.1 * 1/f ;
26.     Nt = floor(tf/dt) ;
27.     for j=1:Nt
28.         t=dt*(j-1);
29.         for i=1:N
30.             if j==1, yp1(i)=0; %eq (7)
31.             elseif j==2
32.                 if i==1
33.                     yp1(1) = EI/2/ro*dt^2/dx^4* ( 2*dx^3*p*cos(w*t)/EI ); %eq (8)
34.                 elseif i==2,
35.                     G(i) = g /(1 + g*( abs( ( y(i+1) - y(i-1) )/(2*a) ) )/tu)^2 ;
36.                     q(i) = kc0*y(i)/( 1 + kc0*y(i)/qcu ) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% within the sc
37.                     yp1(i)= y(i) - q(i)/2/ro*dt^2 -EI/2/ro *dt^2/dx^4*(y(i+2) -4*y(i+1)+7*y(i)-4*y(i-1));
38.                     %eq(23)>> 6
39.                 elseif i==N, yp1(i)=0;
40.                     %eq(25)>>10
41.                 else
42.                     G(i) = g /(1 + g*( abs( ( y(i+1) - y(i-1) )/(2*a) ) )/tu)^2 ;
43.                     if ( ( 0  <= rem( abs(x(i)),SL*x_f ) ) && ( rem(abs(x(i)),SL*x_f ) <= DL/2*x_f  ) ) || ( ( ( SL-DL/2)*x_f  <= rem( abs(x(i)),SL*x_f ) ) && ( rem( abs(x(i)),SL*x_f ) <= SL*x_f ) )
44.                         q(i) = kc0*y(i)/( 1 + kc0*y(i)/qcu ) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% within the sc
45.                     else
46.                         q(i) = ks0*y(i)/( u*(1 + ks0*y(i)/qu )) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% Within the soil
47.                     end
48.
49.                     yp1(i) = y(i) - q(i)/2/ro*dt^2 - EI/2/ro*dt^2/dx^4*(y(i+2)-4*y(i+1)+6*y(i)-4*y(i-1)+y(i-2)) ;
50.                     %eq(20)>>4
51.                 end
52.             else
53.                 if i==1
54.                     G(1) = g ;
55.                     q(1) = kc0*y(1)/( 1 + kc0*y(1)/qcu ) - G(1)*H*( 2*y(2) - 2*y(1) )/a^2 ;% within the sc
56.
57.                     yp1(1) = 2*y(1) - ym1(1) - q(1) /ro*dt^2 - EI/ro * dt^2/dx^4*(2*y(3)-8*y(2) + 6*y(1)-2*dx^3*p/EI*cos(w*t) );
58.                     %eq (7)
59.                 elseif i==2,
60.                     G(i) = g /(1 + g*( abs( ( y(i+1) - y(i-1) )/(2*a) ) )/tu)^2 ;
61.                     q(i) = kc0*y(i)/( 1 + kc0*y(i)/qcu ) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% within the sc
62.
63.                     yp1(i) = 2*y(i) - ym1(i) - q(i)/ro *dt^2 - EI/ro*dt^2/dx^4*(y(i+2) -4*y(i+1)+7*y(i)-4*y(i-1) );
64.                     %eq(22)>>5
65.                 elseif i==N,
66.                     G(i) = g /(1 + g*( abs( ( y(i+1) - y(i-1) )/(2*a) ) )/tu)^2 ;
67.                     if ( ( 0  <= rem( abs(x(i)),SL*x_f ) ) && ( rem(abs(x(i)),SL*x_f ) <= DL/2*x_f  ) ) || ( ( ( SL-DL/2)*x_f  <= rem( abs(x(i)),SL*x_f ) ) && ( rem( abs(x(i)),SL*x_f ) <= SL*x_f ) )
68.                         q(i) = kc0*y(i)/( 1 + kc0*y(i)/qcu ) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% within the sc
69.                     else
70.                         q(i) = ks0*y(i)/( u*(1 + ks0*y(i)/qu )) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% Within the soil
71.                     end
72.                     yp1(i)= 2*y(i) - ym1(i) - q(i)/ro*dt^2 - EI/ro*dt^2/dx^4*(y(i-2) - 4*y(i-1)+5*y(i) );
73.                     %eq(24)>>9
74.                 else
75.                     G(i) = g /(1 + g*( abs( ( y(i+1) - y(i-1) )/(2*a) ) )/tu)^2 ;
76.                     if ( ( 0  <= rem( abs(x(i)),SL*x_f ) ) && ( rem(abs(x(i)),SL*x_f ) <= DL/2*x_f  ) ) || ( ( ( SL-DL/2)*x_f  <= rem( abs(x(i)),SL*x_f ) ) && ( rem( abs(x(i)),SL*x_f ) <= SL*x_f ) )
77.                         q(i) = kc0*y(i)/( 1 + kc0*y(i)/qcu ) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% within the sc
78.                     else
79.                         q(i) = ks0*y(i)/( u*(1 + ks0*y(i)/qu )) - G(i)*H*( y(i-1) - 2*y(i) + y(i+1) )/a^2 ;% Within the soil
80.                     end
81.                     yp1(i)=2*y(i)-ym1(i) - q(i)/ro*dt^2 - EI/ro*dt^2/dx^4*( y(i+2) -4*y(i+1)+6*y(i)-4*y(i-1)+y(i-2));
82.                     %eq(17)>>1
83.                 end
84.             end
85.         end
86.         ym1 = y;
87.         y = yp1;
88.
89.     end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top