Shear and adiabatic heating in the p-nodes.

%This program complished the process of 2D Eulerian advection with method
%of allocating memory temporarily.
clear all;clf;


% Set parameters
Nt = 1;
output_interval = 1;
Lx = 500000; 
Ly = 500000;
Nx = 51; 
Ny = 51; 
Nx1 = Nx+1;
Ny1 = Ny+1;
dx = Lx / (Nx-1);
dy = Ly / (Ny-1);
% Basic E-nodes.
x = 0:dx:Lx;
y = 0:dy:Ly;
% Vx staggered nodes.(Nx_node * Ny_node+1)
xvx = 0:dx:Lx;
yvx = -dy/2 : dy : Ly+dy/2;
% Vy staggered nodes.(Nx_node+1 * Ny_node)
xvy = -dx/2 : dx : Lx+dx/2;
yvy = 0:dy:Ly;
% P staggered nodes.(Nx_node+1 * Ny_node+1)
xp = -dx/2 : dx : Lx+dx/2;
yp = -dy/2 : dy : Ly+dy/2;
vx_matrix = zeros(Ny1,Nx1);
vy_matrix = zeros(Ny1,Nx1);
pvx = zeros(Ny1,Nx1);
pvy = zeros(Ny1,Nx1);
vxp=zeros(Ny1,Nx1); 
vyp=zeros(Ny1,Nx1);
p_matrix = zeros(Ny1,Nx1);
g_y = 10;
ETA = zeros(Ny,Nx);%density
RHO = zeros(Ny,Nx);%viscosity
% Amount of markers.
num_per_cell = 4;
Nx_marker = (Nx-1)*num_per_cell;
Ny_marker = (Ny-1)*num_per_cell;
marknum=Nx_marker*Ny_marker;
% Step-length among L-marker points
x_m_step = Lx/Nx_marker;
y_m_step = Ly/Ny_marker;
% Marker points.
xm=zeros(1,marknum); % Horizontal coordinates, m
ym=zeros(1,marknum);
rhom=zeros(1,marknum); % Density, kg/m^3
etam=zeros(1,marknum);
[X,Y] = meshgrid( x_m_step/2:x_m_step:Lx - x_m_step/2 , y_m_step/2:y_m_step:Ly - y_m_step/2 );
vis_marker = zeros(Ny_marker,Nx_marker);
den_marker = zeros(Ny_marker,Nx_marker);
den_plume = 3200;
radius = 100000;
den_mantle = 3300;
x_mid = Lx/2;
y_mid = Ly/2;
RHOVY=zeros(Ny1,Nx1); % Density in vy-nodes, kg/m^3

x_node = zeros(Ny_marker,Nx_marker);
y_node = zeros(Ny_marker,Nx_marker);

in_node = zeros(Ny+1,Nx+1);
in_vx = zeros(Ny+1,Nx+1);
in_vy = zeros(Ny+1,Nx+1);
in_p = zeros(Ny+1,Nx+1);
vx_markers = zeros(Ny_marker,Nx_marker);
vy_markers = zeros(Ny_marker,Nx_marker);


% Parameters and matrixs of exercise 9.3
srxy = zeros(Ny,Nx);
dsxy = zeros(Ny,Nx);
srxx = zeros(Ny,Nx);
dsxx = zeros(Ny,Nx);
pxyxy = zeros(Ny,Nx);
pxxxx = zeros(Ny,Nx);
Hs = zeros(Ny1,Nx1);
T = 1300;%K
aef = 3e-5;
Ha = zeros(Ny1,Nx1);

m = 1;
for jm=1:1:Nx_marker
    for im=1:1:Ny_marker
        % Define marker coordinates
        xm(m)=x_m_step/2+(jm-1)*x_m_step+(rand-0.5)*x_m_step;
        ym(m)=y_m_step/2+(im-1)*y_m_step+(rand-0.5)*y_m_step;
        % Marker properties
        rmark=((xm(m)-x_mid)^2+(ym(m)-y_mid)^2)^0.5;
        if(rmark<=radius)
            rhom(m)=3200; % Mantle density
            etam(m)=1e+20; % Mantle viscosity
        elseif(ym(m)<=0.2*Ly)
            rhom(m)=1; % Sticky density
            etam(m)=1e+17; % Sticky viscosity
        else
            rhom(m)=3300; % Plume density
            etam(m)=1e+21; % Plume viscosity
        end
        % Update marker counter
        m=m+1;
    end
end

Kcont = 1e+21/dx;

dt=0e+12; % initial timestep

% Unknown_number is (Ny+1)*(Nx+1)*3
unknown_number = (Nx+1)*(Ny+1)*3;
coe = sparse(unknown_number,unknown_number);
R = zeros(unknown_number, 1);
 
% Boundary condition : free slip = -1 ; no slip = 1
bc_left   = -1;
bc_right  = -1;
bc_top    = -1;
bc_bottom = -1;
dxymax = 0.5;
%  Process of 2D transport
figure(1);
for t = 1:Nt  
    w_m_x_node = zeros(Ny,Nx);
    den_w_m = zeros(Ny,Nx);
    vis_w_m = zeros(Ny,Nx);
    RHOYSUM=zeros(Ny1,Nx1);
    WTYSUM=zeros(Ny1,Nx1);
    % "drunken sailor"  instability appear if dt = 0 at the next line.
%     dt = 0;
    % Interpolate density and viscosity from markers to BASIC NODES.
    % Bisection algorithm
    % Searching the coordinate of upper-left node of every L-marker.
    for m=1:1:marknum
        % Define i,j indexes for the upper left node
        j=fix((xm(m)-x(1))/dx)+1;
        i=fix((ym(m)-y(1))/dy)+1;
        if(j<1)
            j=1;
        elseif(j>Nx-1)
            j=Nx-1;
        end
        if(i<1)
            i=1;
        elseif(i>Ny-1)
            i=Ny-1;
        end
        % Compute distances
        dis_x = xm(m)-x(j);
        dis_y = ym(m)-y(i);
        wtmij = (1-dis_x/dx)*(1-dis_y/dy);
        wtmi1j = (dis_x/dx)*(1-dis_y/dy);
        wtmij1 = (1-dis_x/dx)*(dis_y/dy);
        wtmi1j1 = (dis_x/dx)*(dis_y/dy);
        % |       |
        %     * * * * *   Discard the m-point on the right.
        w_m_x_node(i,j) = w_m_x_node(i,j) + wtmij;
        den_w_m(i,j) = den_w_m(i,j) + rhom(m).*wtmij;
        vis_w_m(i,j) = vis_w_m(i,j) + etam(m).*wtmij;
        
        w_m_x_node(i+1,j) = w_m_x_node(i+1,j) + wtmi1j;
        den_w_m(i+1,j) = den_w_m(i+1,j) + rhom(m).*wtmi1j;
        vis_w_m(i+1,j) = vis_w_m(i+1,j) + etam(m).*wtmi1j;
        
        w_m_x_node(i,j+1) = w_m_x_node(i,j+1) + wtmij1;
        den_w_m(i,j+1) = den_w_m(i,j+1) + rhom(m).*wtmij1;
        vis_w_m(i,j+1) = vis_w_m(i,j+1) + etam(m).*wtmij1;
        
        w_m_x_node(i+1,j+1) = w_m_x_node(i+1,j+1) + wtmi1j1;
        den_w_m(i+1,j+1) = den_w_m(i+1,j+1) + rhom(m).*wtmi1j1;
        vis_w_m(i+1,j+1) = vis_w_m(i+1,j+1) + etam(m).*wtmi1j1;
            
        % Density interpolation to vy-nodes
        % Define i,j indexes for the upper left node
        % Vy staggered nodes.(Nx_node+1 * Ny_node)
        j=fix((xm(m)-xvy(1))/dx)+1;
        i=fix((ym(m)-yvy(1))/dy)+1;
        if(j<1)
            j=1;
        elseif(j>Nx)
            j=Nx;
        end
        if(i<1)
            i=1;
        elseif(i>Ny-1)
            i=Ny-1;
        end
        % Compute distances
        dxmj=xm(m)-xvy(j);
        dymi=ym(m)-yvy(i);
        % Compute weights
        wtmij=(1-dxmj/dx)*(1-dymi/dy);
        wtmi1j=(1-dxmj/dx)*(dymi/dy);
        wtmij1=(dxmj/dx)*(1-dymi/dy);
        wtmi1j1=(dxmj/dx)*(dymi/dy);
        % Update properties
        % i,j Node
        RHOYSUM(i,j)=RHOYSUM(i,j)+rhom(m)*wtmij;
        WTYSUM(i,j)=WTYSUM(i,j)+wtmij;
        % i+1,j Node
        RHOYSUM(i+1,j)=RHOYSUM(i+1,j)+rhom(m)*wtmi1j;
        WTYSUM(i+1,j)=WTYSUM(i+1,j)+wtmi1j;
        % i,j+1 Node
        RHOYSUM(i,j+1)=RHOYSUM(i,j+1)+rhom(m)*wtmij1;
        WTYSUM(i,j+1)=WTYSUM(i,j+1)+wtmij1;
        % i+1,j+1 Node
        RHOYSUM(i+1,j+1)=RHOYSUM(i+1,j+1)+rhom(m)*wtmi1j1;
        WTYSUM(i+1,j+1)=WTYSUM(i+1,j+1)+wtmi1j1;
    end
    % Computing the value of density and viscosity of E-nodes. 
    for j = 1:Nx
        for i = 1:Ny
            if(w_m_x_node(i,j)>0)
                % Not need RHO on the basic nodes.Just Y nodes.
                RHO(i,j) = den_w_m(i,j)/w_m_x_node(i,j);
                ETA(i,j) = vis_w_m(i,j)/w_m_x_node(i,j);
            end
        end
    end
    for j=1:1:Nx1
        for i=1:1:Ny1
            if(WTYSUM(i,j)>0)
                RHOVY(i,j)=RHOYSUM(i,j)/WTYSUM(i,j);
            end
        end
    end
    ETAP = zeros(Ny+1,Nx+1);
    % Computing the harmonic density.
    for j = 2:Nx
        for i = 2:Ny
            % Range of var_vis_n is 2:Ny 2:Nx.
            ETAP(i,j) = 1/( (1/ETA(i,j) + 1/ETA(i-1,j) + 1/ETA(i,j-1) + 1/ETA(i-1,j-1) )/4 );
        end
    end
    for j=1:1:Nx1
        for i=1:1:Ny1
            % Define global indexes in algebraic space
            kvx=((j-1)*Ny1+i-1)*3+1; % Vx
            kvy=kvx+1; % Vy
            kpm=kvx+2; % P
            
            % Vx equation External points
            if(i==1 || i==Ny1 || j==1 || j==Nx || j==Nx1)
                % Boundary Condition
                % 1*Vx=0
                coe(kvx,kvx)=1; % Left part
                R(kvx)=0; % Right part
                % Top boundary
                if(i==1 && j>1 && j<Nx)
                    coe(kvx,kvx+3)=bc_top; % Left part
                end
                % Bottom boundary
                if(i==Ny1 && j>1 && j<Nx)
                    coe(kvx,kvx-3)=bc_bottom; % Left part
                end
            else
                % Internal points: x-Stokes eq.
                % ETA*(d2Vx/dx^2+d2Vx/dy^2)-dP/dx=0
                %            Vx2
                %             |
                %        Vy1  |  Vy3
                %             |
                %     Vx1-P1-Vx3-P2-Vx5
                %             |
                %        Vy2  |  Vy4
                %             |
                %            Vx4
                %
                % Viscosity points
                ETA1=ETA(i-1,j);
                ETA2=ETA(i,j);
                ETAP1=ETAP(i,j);
                ETAP2=ETAP(i,j+1);
                % Left part
                coe(kvx,kvx-Ny1*3)=2*ETAP1/dx^2; % Vx1
                coe(kvx,kvx-3)=ETA1/dy^2; % Vx2
                coe(kvx,kvx)=-2*(ETAP1+ETAP2)/dx^2-(ETA1+ETA2)/dy^2; % Vx3
                coe(kvx,kvx+3)=ETA2/dy^2; % Vx4
                coe(kvx,kvx+Ny1*3)=2*ETAP2/dx^2; % Vx5
                coe(kvx,kvy)=-ETA2/dx/dy;  % Vy2
                coe(kvx,kvy+Ny1*3)=ETA2/dx/dy;  % Vy4
                coe(kvx,kvy-3)=ETA1/dx/dy;  % Vy1
                coe(kvx,kvy+Ny1*3-3)=-ETA1/dx/dy;  % Vy3
                coe(kvx,kpm)=Kcont/dx; % P1
                coe(kvx,kpm+Ny1*3)=-Kcont/dx; % P2
                % Right part
                R(kvx)=0;
            end
            
            % Vy equation External points
            if(j==1 || j==Nx1 || i==1 || i==Ny || i==Ny1)
                % Boundary Condition
                % 1*Vy=0
                coe(kvy,kvy)=1; % Left part
                R(kvy)=0; % Right part
                % Left boundary
                if(j==1 && i>1 && i<Ny)
                    coe(kvy,kvy+3*Ny1)=bc_left; % Left part
                end
                % Right boundary
                if(j==Nx1 && i>1 && i<Ny)
                    coe(kvy,kvy-3*Ny1)=bc_right; % Left part
                end
            else
                % Internal points: y-Stokes eq.
                % ETA*(d2Vy/dx^2+d2Vy/dy^2)-dP/dy=-RHO*gy
                %            Vy2
                %             |
                %         Vx1 P1 Vx3
                %             |
                %     Vy1----Vy3----Vy5
                %             |
                %         Vx2 P2 Vx4
                %             |
                %            Vy4
                %
                % Viscosity points
                % Viscosity points
                ETA1=ETA(i,j-1);
                ETA2=ETA(i,j);
                ETAP1=ETAP(i,j);
                ETAP2=ETAP(i+1,j);
                dRHOdx=(RHOVY(i,j+1)-RHOVY(i,j-1))/2/dx;
                dRHOdy=(RHOVY(i+1,j)-RHOVY(i-1,j))/2/dy;

                % Left part
                coe(kvy,kvy-Ny1*3)=ETA1/dx^2; % Vy1
                coe(kvy,kvy-3)=2*ETAP1/dy^2; % Vy2
                coe(kvy,kvy)=-2*(ETAP1+ETAP2)/dy^2-...
                    (ETA1+ETA2)/dx^2 - dRHOdy*g_y*dt; % Vy3
                coe(kvy,kvy+3)=2*ETAP2/dy^2; % Vy4
                coe(kvy,kvy+Ny1*3)=ETA2/dx^2; % Vy5
                coe(kvy,kvx)=-ETA2/dx/dy -dRHOdx*g_y*dt/4; %Vx3
                coe(kvy,kvx+3)=ETA2/dx/dy -dRHOdx*g_y*dt/4; %Vx4
                coe(kvy,kvx-Ny1*3)=ETA1/dx/dy -dRHOdx*g_y*dt/4; %Vx1
                coe(kvy,kvx+3-Ny1*3)=-ETA1/dx/dy -dRHOdx*g_y*dt/4; %Vx2
                coe(kvy,kpm)=Kcont/dy; % P1
                coe(kvy,kpm+3)=-Kcont/dy; % P2
                
                % Right part
                R(kvy)=-RHOVY(i,j)*g_y;
            end
            
            % P equation External points
            if(i==1 || j==1 || i==Ny1 || j==Nx1 ||...
                    (i==2 && j==2))
                % Boundary Condition
                % 1*P=0
                coe(kpm,kpm)=1; % Left part
                R(kpm)=0; % Right part
                % Real BC
                if(i==2 && j==2)
                    coe(kpm,kpm)=1*Kcont; %Left part
                    R(kpm)=1e+9; % Right part
                end
            else
                % Internal points: continuity eq.
                % dVx/dx+dVy/dy=0
                %            Vy1
                %             |
                %        Vx1--P--Vx2
                %             |
                %            Vy2
                %
                % Left part
                coe(kpm,kvx-Ny1*3)=-1/dx; % Vx1
                coe(kpm,kvx)=1/dx; % Vx2
                coe(kpm,kvy-3)=-1/dy; % Vy1
                coe(kpm,kvy)=1/dy; % Vy2
                % Right part
                R(kpm)=0;
            end
            
        end
    end
    % Computing the solution U:include vx,vy and p
     U = coe \ R;
     
    for j=1:1:Nx1
        for i=1:1:Ny1
            % Define global indexes in algebraic space
            kvx=((j-1)*Ny1+i-1)*3+1; % Vx
            kvy=kvx+1; % Vy
            kpm=kvx+2; % P
            % Reload solution
            vx_matrix(i,j)=U(kvx);
            vy_matrix(i,j)=U(kvy);
            p_matrix(i,j)=U(kpm)*Kcont;
        end
    end
    % Computing strain rate and deviatoric stress components.
    for j = 1:Nx
        for i = 1:Ny
            srxy(i,j) = (1/2)*( (vx_matrix(i+1,j)-vx_matrix(i,j))/dy+(vy_matrix(i,j+1)-vy_matrix(i,j))/dx ) ;
            dsxy(i,j) = 2*ETA(i,j)*srxy(i,j);
            if(j~=1)
                srxx(i,j) = ( vx_matrix(i,j)-vx_matrix(i,j-1) )/dx;
            end
            dsxx(i,j) = 2*ETAP(i,j)*srxx(i,j);
            if(i~=1 && j~=1)
                kxyxyij = srxy(i,j)*dsxy(i,j);
                kxyxyi1j = srxy(i-1,j)*dsxy(i-1,j);
                kxyxyij1 = srxy(i,j-1)*dsxy(i,j-1);
                kxyxyi1j1 = srxy(i-1,j-1)*dsxy(i-1,j-1);
                pxyxy(i,j) = (1/4)*(kxyxyij+kxyxyi1j+kxyxyij1+kxyxyi1j1);
                kxxxxij = srxx(i,j)*dsxx(i,j);
                kxxxxi1j = srxx(i-1,j)*dsxx(i-1,j);
                kxxxxij1 = srxx(i,j-1)*dsxx(i,j-1);
                kxxxxi1j1 = srxx(i-1,j-1)*dsxx(i-1,j-1);
                pxxxx(i,j) = (1/4)*(kxxxxij+kxxxxi1j+kxxxxij1+kxxxxi1j1);
            end
            % Shear heating
            Hs(i,j) = 2*( pxyxy(i,j) + pxxxx(i,j) ) ;
        end
    end
    % Adiabatic heating
    for j = 1:Nx1
        for i = 2:Ny1
            kvy = ( vy_matrix(i,j)+vy_matrix(i-1,j) )/2;
            kRHO = ( RHOVY(i,j)+RHOVY(i-1,j) )/2;
            Ha(i,j) = T*aef*g_y*kvy*kRHO;
        end
    end
    
    
    % Computing average velocity components at cell centres.
    for j = 1:Nx1
        for i = 1:Ny1
            %  Computing internal Vx and Vy at the P nodes.
            if (j~=1 && i ~=1 && j~=Nx1 && i~=Ny1)
                % Vx at internal nodes.
                pvx(i,j) = ( vx_matrix(i,j) + vx_matrix(i,j-1) )./2;           
                % Vy at internal nodes.
                pvy(i,j) = ( vy_matrix(i,j) + vy_matrix(i-1,j) )./2;           
            end
        end
    end
    % Why it is free-slip ??
    % Applying free-slip boundary condition for velocity in the P nodes .
    for j = 1:Nx1
        for i = 1:Ny1
            %    j=1 j=2 j=3 ... j=Nx_node+1
            % i=1
            % i=2
            % i=3
            % ...
            % i=Ny_node+1
            if ( i==1 || i==Ny )
                % Top and Bottom
                % Vx
                if(j>1 && j<Nx1)
                    pvx(i,j) = pvx(i+1,j);
                end
                % Vy
                pvy(i,j) = pvy(i+1,j);
            end
            if (j==1 || j==Nx)
                % Right and Left
                % Vy
                if(i>1 && i< Ny1)
                    pvy(i,j) = pvy(i,j+1);
                end
                % Vx
                pvx(i,j) = pvx(i,j+1);
            end
        end
    end
    % Define time steps.
    dt = 1e+36;
    maxvx=max(max(abs(vx_matrix)));
    maxvy=max(max(abs(vy_matrix)));
    if(dt*maxvx>dxymax*dx)
        dt=dxymax*dx/maxvx;
    end
    if(dt*maxvy>dxymax*dy)
        dt=dxymax*dy/maxvy;
    end
    % The classical fourth-order Runge-Kutta scheme.
    % Define Vxa(Vya)1,Vxb(Vyb)2,Vxc(Vyc)3,Vxd(Vyd)4.
    vxm=zeros(4,1);
    vym=zeros(4,1);
    % Interpolate vx and vy components from E-nodes to L-markers.
    for m=1:1:marknum
        % Reserve original coordinate.
        xa = xm(m);
        ya = ym(m);
        for rk = 1:1:4
            % Interpolate vx from P-staggered nodes to L-markers.
            % P staggered nodes.(Nx_node+1 * Ny_node+1)
            j=fix((xm(m)-xp(1))/dx)+1;
            i=fix((ym(m)-yp(1))/dy)+1;
            if(j<1)
                j=1;
            elseif(j>Nx-1)
                j=Nx-1;
            end
            % Vx staggered nodes.(Nx_node * Ny_node+1)
            if(i<1)
                i=1;
            elseif(i>Ny)
                i=Ny;
            end
            % Compute distances
            dis_x=xm(m)-xp(j);
            dis_y=ym(m)-yp(i);
            % Compute weights
            wtmij=(1-dis_x/dx)*(1-dis_y/dy);
            wtmi1j=(dis_x/dx)*(1-dis_y/dy);
            wtmij1=(1-dis_x/dx)*(dis_y/dy);
            wtmi1j1=(dis_x/dx)*(dis_y/dy);
            % Compute vx velocity
            vxm(rk) = pvx(i,j)*wtmij+pvx(i+1,j)*wtmi1j...
                     +pvx(i,j+1)*wtmij1+pvx(i+1,j+1)*wtmi1j1;
            vym(rk) = pvy(i,j)*wtmij+pvy(i+1,j)*wtmi1j...
                     +pvy(i,j+1)*wtmij1+pvy(i+1,j+1)*wtmi1j1;
            % Interpolate vx from Vx-staggered nodes to L-markers.
            % Define i,j indexes for the upper left node
            j=fix((xm(m)-xvx(1))/dx)+1;
            i=fix((ym(m)-yvx(1))/dy)+1;
            if(j<1)
                j=1;
            elseif(j>Nx-1)
                j=Nx-1;
            end
            % Vx staggered nodes.(Nx_node * Ny_node+1)
            if(i<1)
                i=1;
            elseif(i>Ny)
                i=Ny;
            end
            % Compute distances
            dis_x=xm(m)-xvx(j);
            dis_y=ym(m)-yvx(i);
            % Compute weights
            wtmij=(1-dis_x/dx)*(1-dis_y/dy);
            wtmi1j=(dis_x/dx)*(1-dis_y/dy);
            wtmij1=(1-dis_x/dx)*(dis_y/dy);
            wtmi1j1=(dis_x/dx)*(dis_y/dy);
            % Compute vx velocity
            vxcij = (1/3)*(2*vx_matrix(i,j));
            vxci1j = (1/3)*(2*vx_matrix(i+1,j));
            vxcij1 = (1/3)*(2*vx_matrix(i,j+1));
            vxci1j1 = (1/3)*(2*vx_matrix(i+1,j+1));
            vxm(rk) = (1/3)*vxm(rk) + (vxcij*wtmij+vxci1j*wtmi1j...
                     +vxcij1*wtmij1+vxci1j1*wtmi1j1);
            
            % Interpolate vy
            % Define i,j indexes for the upper left node
            j=fix((xm(m)-xvy(1))/dx)+1;
            i=fix((ym(m)-yvy(1))/dy)+1;
            if(j<1)
                j=1;
            elseif(j>Nx)
                j=Nx;
            end
            if(i<1)
                i=1;
            elseif(i>Ny-1)
                i=Ny-1;
            end
            % Compute distances
            dis_x=xm(m)-xvy(j);
            dis_y=ym(m)-yvy(i);
            % Compute weights
            wtmij=(1-dis_x/dx)*(1-dis_y/dy);
            wtmi1j=(1-dis_x/dx)*(dis_y/dy);
            wtmij1=(dis_x/dx)*(1-dis_y/dy);
            wtmi1j1=(dis_x/dx)*(dis_y/dy);
            % Compute vx velocity
            vycij = (1/3)*(2*vy_matrix(i,j));
            vyci1j = (1/3)*(2*vy_matrix(i+1,j));
            vycij1 = (1/3)*(2*vy_matrix(i,j+1));
            vyci1j1 = (1/3)*(2*vy_matrix(i+1,j+1));
            vym(rk) = (1/3)*vym(rk) + (vycij*wtmij+vyci1j*wtmi1j...
                     +vycij1*wtmij1+vyci1j1*wtmi1j1);
                 
            if(rk == 1 || rk == 2)
                xm(m) = vx_matrix(rk)*dt/2 + xa;
                ym(m) = vy_matrix(rk)*dt/2 + ya;
            elseif(rk == 3)
                xm(m) = vx_matrix(rk)*dt + xa;
                ym(m) = vy_matrix(rk)*dt + ya;
            end
        end
        % Return original coordinate.
        xm(m) = xa;
        ym(m) = ya;
        % Move markers
        vxefft = (1/6)*(vxm(1)+2*vxm(2)+2*vxm(3)+vxm(4));
        vyefft = (1/6)*(vym(1)+2*vym(2)+2*vym(3)+vym(4));        
        xm(m)=xm(m)+dt*vxefft;
        ym(m)=ym(m)+dt*vyefft;
    end
    figure(1);
    colormap('Jet');clf
    pcolor(x,y,log10(ETA));
    hold on;
    %P nodes:(Nx+1 * Ny+1).
    quiver(xp(1:2:Nx1),yp(1:2:Ny1),pvx(1:2:Ny1,1:2:Nx1),pvy(1:2:Ny1,1:2:Nx1),'w')
    xlabel('Horizontal(m)')
    ylabel('Vertical(m)')
    caxis([17 21]);
    set(gca,'xaxislocation','top');
    set (gca,'YDir','reverse')
    shading interp;
    colorbar
    title(['log10(ETA) after ',num2str(t-1),' deltat'])
    pause(0.01);
    % Strain rate component:srxy
    figure(2);
    colormap('Jet');
    pcolor(x(1:Ny),y(1:Nx),srxy);
    xlabel('Horizontal(m)')
    ylabel('Vertical(m)')
    set(gca,'xaxislocation','top');
    set (gca,'YDir','reverse')
    shading interp;
    colorbar
    title(['srxy after ',num2str(t-1),' deltat'])
    pause(0.01);
    % Deviatoric stress component:srxy
    figure(3);
    colormap('Jet');
    pcolor(x(1:Ny),y(1:Nx),dsxy);
    xlabel('Horizontal(m)')
    ylabel('Vertical(m)')
    set(gca,'xaxislocation','top');
    set (gca,'YDir','reverse')
    shading interp;
    colorbar
    title(['dsxy after ',num2str(t-1),' deltat'])
    pause(0.01);
    % Deviatoric stress component:srxy
    figure(4);
    colormap('Jet');
    pcolor(x(1:Ny),y(1:Nx),srxx);
    xlabel('Horizontal(m)')
    ylabel('Vertical(m)')
    set(gca,'xaxislocation','top');
    set (gca,'YDir','reverse')
    shading interp;
    colorbar
    title(['srxx after ',num2str(t-1),' deltat'])
    pause(0.01);
    % Deviatoric stress component:srxy
    figure(5);
    colormap('Jet');
    pcolor(x(1:Ny),y(1:Nx),dsxx);
    xlabel('Horizontal(m)')
    ylabel('Vertical(m)')
    set(gca,'xaxislocation','top');
    set (gca,'YDir','reverse')
    shading interp;
    colorbar
    title(['dsxx after ',num2str(t-1),' deltat'])
    pause(0.01);  
    % Hs
    figure(6);
    colormap('Jet');
    pcolor(xp,yp,Hs);
    xlabel('Horizontal(m)')
    ylabel('Vertical(m)')
    set(gca,'xaxislocation','top');
    set (gca,'YDir','reverse')
    shading interp;
    colorbar
    title(['Hs after ',num2str(t-1),' deltat'])
    pause(0.01);  
    % Hs
    figure(7);
    colormap('Jet');
    pcolor(xp,yp,Ha);
    xlabel('Horizontal(m)')
    ylabel('Vertical(m)')
    set(gca,'xaxislocation','top');
    set (gca,'YDir','reverse')
    shading interp;
    colorbar
    title(['Ha after ',num2str(t-1),' deltat'])
    pause(0.01);  
        
end

 

 

 

 

 

 

 

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
To find the tensile capacity P1 for the single-shear lap connection, we need to calculate the strength of the bolts and the strength of the plate. Strength of Bolts: The ultimate strength of the bolts is determined by the minimum of the following: 1. Bolt shear strength 2. Bolt bearing strength 3. Bolt tensile strength Bolt Shear Strength: The bolt shear strength is given by: Vn = 0.6*Fy*Ab where Vn is the nominal shear strength of the bolt, Fy is the yield strength of the bolt material, and Ab is the cross-sectional area of the bolt. Substituting the given values, we get: Vn = 0.6*90*0.44 = 23.76 kips Bolt Bearing Strength: The bolt bearing strength is given by: Rn = 2.4*dt*t*Fyb where Rn is the nominal bearing strength of the bolt, dt is the thickness of the connected material, t is the thickness of the bolt head or nut, and Fyb is the yield strength of the bolt material. Substituting the given values, we get: Rn = 2.4*0.75*0.75*90 = 121.5 kips Bolt Tensile Strength: The bolt tensile strength is given by: Tn = 0.75*Fub*Ab where Tn is the nominal tensile strength of the bolt, Fub is the ultimate strength of the bolt material, and Ab is the cross-sectional area of the bolt. Substituting the given values, we get: Tn = 0.75*120*0.44 = 29.7 kips The minimum of the above three values is the ultimate strength of the bolt, which is 23.76 kips. Strength of Plate: The ultimate strength of the plate is determined by the minimum of the following: 1. Plate tensile strength 2. Plate bearing strength 3. Plate yielding strength Plate Tensile Strength: The plate tensile strength is given by: Pu = Fu*Ap where Pu is the ultimate strength of the plate, Fu is the ultimate strength of the plate material, and Ap is the effective net area of the plate. The effective net area of the plate is given by: Ap = An - n*d*(dh + 0.5*t) where An is the gross area of the plate, n is the number of bolts, d is the bolt diameter, dh is the diameter of the hole, and t is the thickness of the connected material. Substituting the given values, we get: An = 12*1 = 12 in² n = 4 d = 0.75 in dh = 0.8125 in (assuming 1/16 in oversize hole) t = 1 in Ap = 12 - 4*0.75*(0.8125 + 0.5*1) = 5.25 in² Substituting the given values, we get: Pu = 58*5.25 = 304.5 kips Plate Bearing Strength: The plate bearing strength is given by: Rn = 2.4*dt*t*Fyb where Rn is the nominal bearing strength of the plate, dt is the thickness of the connected material, t is the thickness of the bolt head or nut, and Fyb is the yield strength of the bolt material. Substituting the given values, we get: Rn = 2.4*1*0.75*36 = 64.8 kips Plate Yielding Strength: The plate yielding strength is given by: Py = 0.9*Fy*Ag where Py is the yielding strength of the plate, Fy is the yield strength of the plate material, and Ag is the gross area of the plate. Substituting the given values, we get: Ag = 12*1 = 12 in² Py = 0.9*36*12 = 388.8 kips The minimum of the above three values is the ultimate strength of the plate, which is 64.8 kips. Now, we can find the tensile capacity P1 of the connection by taking the minimum of the ultimate strength of the bolts and the ultimate strength of the plate. P1 = min(23.76, 64.8) = 23.76 kips Therefore, the tensile capacity P1 for the single-shear lap connection is 23.76 kips.

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值