3.1 Preparación del modelo con GEFDyn

Existen varias opciones para modelar el problema dependiendo del esquema adoptado para inicializar los esfuerzos efectivos al interior del suelo. Ya que se desea estudiar el comportamiento inelástico del material, la etapa de inicialización es fundamental:

En términos de las opciones de GEFDyn, la primera opción consiste en efectuar un primer cálculo estático empleando MODEX=1 tomando en cuenta las fuerzas de gravedad, seguido de un cálculo dinámico con MODEX=6, es decir, retomando el cálculo estático poniendo en cero deformaciones y desplazamientos.

La segunda opción sería efectuar un único cálculo dinámico con MODEX=1, definiendo un ``buen'' estado inicial de esfuerzos por ejemplo empleando la opción NCOUCH. Si se emplea la opción (IDESEQ=1) el programa asumirá que el estado proporcionado es correcto y guardará las fuerzas internas asociadas para ``equilibrarlas'' (agregarlas al ensamble del vector de fuerzas nodales en cada paso de carga). Luego, se puede efectuar directamente el análisis dinámico sin gravedad.

La primera opción tiene la desventaja que requiere dos cálculos, pero es más versátil en el sentido que durante la fase estática se pueden modelar cualquier tipo de carga: cargas inducidas, excavaciones, rellenos, variaciones de napas freáticas, etc. La segunda opción tiene la ventaja de requerir un único cálculo, pero está limitada a la posibilidad de definir un estado inicial de tensiones efectivas razonable. Ya que la primera opción es más versátil y es la que se emplea mayormente en la práctica, describiremos esta opción en lo que sigue.

Para la construcción del mallado emplearemos una vez más las herramientas de SDTools. Ya que la malla emplear es exactamente la misma (a excepción del elemento paraxial que aparece sólo en la fase dinámica), emplearemos un mismo script variando sólo las condiciones de borde para el caso estático y el dinámico:

clear all
close all
 
%-------------------------------------------------------------------------
%Mesh construction
filename = ('col_ini_2d');

cas=0; %=0 static, =1 dynamic
FEnode=[];  FEelt=[];
%----------------------------------------------------------------------Soil
FEel0=[]; 
FEnode=[1 0 0 0 0. 0 .0
        2 0 0 0 0. 1. .0];
%Layer 1
FEel0=[]; NodeID=[];
NodeID=femesh(['findnode z==0.']);
[Y,I]=sort(FEnode(NodeID,6),'ascend'); NodeID=NodeID(I);
femesh('ObjectBeamLine',NodeID)
femesh('extrude 1 0 0 -30')
femesh('divide 30')
femesh('addsel')
%Bedrock
FEel0=[]; NodeID=[];
NodeID=femesh('findnode z==-30.');    
[Y,I]=sort(FEnode(NodeID,6),'ascend');  NodeID=NodeID(I);   
femesh('ObjectBeamLine',NodeID)
femesh('extrude 1 0 0 -5')
femesh('addsel')
%------------------------------------------------------------------Paraxial
FEel0=[]; 
NodeID=femesh('findnode z==-35.');    
[Y,I]=sort(FEnode(NodeID,6),'ascend');  NodeID=NodeID(I);   
femesh('ObjectBeamLine',NodeID)
femesh('addsel')
%---------------------------------------------------Numeration optimisation 
[FEnode,FEelt]=feutil('optimnodenum',FEnode,FEelt);
model.Node = FEnode;
model.Elt  = FEelt;
model.Elt  = feutil('orient',model);
model.DOF  = feutil('getdof',model);
model=feutil('renumber',model) 
FEnode(:,:)=model.Node;
FEelt(:,:)=model.Elt;
%--------------------------------------------------------Writing file .geom
filename2=[filename,'.geom'];
fid_geom=fopen(filename2,'w');
numnp=length(FEnode(:,1));
%Title writing
fprintf(fid_geom,'********* TEST COLUMN 2D **********\n');
%---------------------------------------------------Truncate to 6th decimal
for inod = 1:numnp
    FEnode(inod,5:7)= round(FEnode(inod,5:7)*1e6)*1e-6;
end
%-------------------------------------------------------------Nodes writing     
for inod = 1:numnp   
    inode=find(FEnode(:,1)==inod);
    ymin=min(FEnode(:,6)); ymax=max(FEnode(:,6)); 
    zmin=min(FEnode(:,7)); zmax=max(FEnode(:,7));
    % id = [x y z rx ry rz p 0 neq]
    id=[1 0 0 1 1 1 1 0 0];     %dry case
    %Equivalent nodes
    if FEnode(inode,6) == ymin 
        NodeID(1)=femesh('findnode x== & y== & z==',...
            FEnode(inode,5),ymax,FEnode(inode,7));
        id(9)=NodeID(1);
    elseif FEnode(inode,6) == ymax
        NodeID(1)=femesh('findnode x== & y== & z==',...
            FEnode(inode,5),ymin,FEnode(inode,7));
        id(9)=NodeID(1);
    end
    %---------------------------------------------------Boundary conditions
    if cas==0 %static
        if ( abs(FEnode(inode,6)-ymin)<1e-5 || abs(FEnode(inode,6)-ymax)<1e-5 ) && ...
                (FEnode(inode,7)>zmin && FEnode(inode,7)<=zmax) ;
            id(2)=1; %Lateral limits(y)
        elseif ( abs(FEnode(inode,6)-ymin)<1e-5 || abs(FEnode(inode,6)-ymax)<1e-5 ) && ...
                abs(FEnode(inode,7)-zmin)<1e-5 
            id(2:3)=[1 1]; %corners
        elseif (FEnode(inode,6)>ymin && FEnode(inode,6)<ymax) && abs(FEnode(inode,7)-zmin)<1e-5 
            id(3)=1; %bottom
        end
    elseif cas==1 %dynamic
         if abs(FEnode(inode,7)-zmin)<1e-5
             id(2:3)=[0 0]; %bottom
         end
    end
    %Writing
    fprintf(fid_geom,'%5i%2i%2i%2i%2i%2i%2i%2i      %10.4f%10.4f%10.4f%5i%6i\n',...
    inod,id(1:7),FEnode(inode,5),FEnode(inode,6),FEnode(inode,7),id(8),id(9));   
end
%-----------------------------------------------------------------Volume 2D
%Parameters for GEFDYN
ITELT=1; IEL=4; IPS=1; EPAI=1.0; MTYP=1; KG=0; 
N_group=[30 1]; %number of elements by group
Group_num=[1 2];
%Quad4
for n_group=1:length(N_group)
    group=Group_num(n_group);
    [ind,ELT] = femesh(['findelt Group ',int2str(group)]);
    for M=1:N_group(n_group)
        fprintf(fid_geom,'%5i%1i%2i%2i%20f%5i%5i%5i%5i%5i%5i\n',...
            M,ITELT,IEL,IPS,EPAI,MTYP,KG,ELT(M+1,2),ELT(M+1,1),ELT(M+1,4),ELT(M+1,3));
    end
end
%------------------------------------------------------------------Paraxial
if cas==1; %dynamic
     N_group=1;     group=3;
     %Parameters
     IEL=2; MTYP=1; KG=0; EPAI=1.0;
     [ind,ELT] = femesh(['findelt Group ',int2str(group)]);
     for M=1:N_group
         fprintf(fid_geom,'%5i%3i%27i%5i%5i%5i%5i%10f\n',...
             M,IEL,MTYP,KG,ELT(M+1,1),ELT(M+1,2),0,EPAI);
     end
end
fprintf(fid_geom,'12345678901234567890123456789012345678901234567890123456789012345678901234567890\n');
%close files
fclose('all')

De acuerdo a la geometría del problema, en el caso estático (cas=0) los bordes laterales sólo pueden desplazarse verticalmente (id(2)=1), mientras que el borde inferior debe estar bloqueado tanto horizontal como verticalmente (id(2:3)=[1 1] ). En este caso hemos supuesto que elementos de 1 [m] de espesor son suficientemente finos para cubrir la gama de frecuencias deseadas. Como la roca posee mayor velocidad de propagación de onda de corte, puede ser mallada más gruesa manteniendo el mismo nivel de precisión. Chequearemos estos valores más abajo, una vez obtenido el perfil de velocidades elástico.

En el caso dinámico, impondremos un cinemática tipo ``viga de corte'' para garantizar la propagación unidimensional. Esto se consigue atando los bordes laterales del mallado mediante la opción de nodos equivalentes. Dos nodos equivalentes poseen los mismos desplazamientos, lo que garantiza la condición buscada. No es razonable emplear elementos laterales absorbentes por dos razones. En primer lugar, se busca representar el comportamiento no lineal del suelo de forma que las condiciones locales elásticas requeridas no se cumplen. En segundo lugar, un elemento absorbente puro no será capaz de discriminar si el campo es incidente o reflejado y tendería a amortiguar el movimiento incluso antes que alcance la superficie y se refleje. En efecto, una segunda alternativa correcta de modelación sería extender el mallado lateralmente hasta una zona razonablemente elástica e imponer sobre esos bordes elementos absorbentes e incidentes sincronizados con los elementos del borde inferior. Finalmente, en el borde inferior se impone un elemento paraxial para reemplazar el semi-espacio infinito1

La malla de elementos finitos generada se ilustra en la Fig.6. La malla cuenta de 31 elementos sólidos 2D (30 del suelo y 1 de la roca), 1 elemento unidimensional (para el paraxial) y 64 nodos. El elemento unidimensional permite automatizar la generación de los elementos paraxiales.

Figura 6: Malla empleada
\begin{figure}\centering
\epsfig{file=figures/tut4_malla_inci.eps,clip=,width=.8\textwidth}
\end{figure}


Notas al pie

... infinito1
notar que en este borde no se impone ninguna condición de borde particular ya que la fuerza de gravedad ya estará equilibrada en el cálculo estático, por lo tanto el modelo no requiere sujeción vertical.


Subsecciones
Esteban Saez 2010-11-24