daily pastebin goal
9%
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
Top