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. OK, I Understand
 
Top