2.1 Preparación del modelo con GEFDyn

Comenzaremos construyendo un malla para el problema. Para efectos de construcción del modelo, emplearemos un ancho unitario y supondremos deformaciones planas para trabajar en dos dimensiones. Para generar la malla con las herramientas de SDT desde Matlab, podemos emplear el script siguiente:

clear all
close all

filename='tut4_parax';%filename
FEnode=[];  FEelt=[]; %reseting default variables
%---------------------------------------------------------Mesh construction
%Solid
FEel0=[]; 
FEnode=[1 0 0 0 0.  0 .0
        2 0 0 0 0. 1. .0];
NodeID=femesh('findnode z==0.0');
[Y,I]=sort(FEnode(NodeID,6),'ascend');  NodeID=NodeID(I);
femesh('ObjectBeamLine',NodeID)
femesh('extrude 1 0 0 -5')
femesh('divide 20 1')
femesh('addsel')
%Parax
NodeID=femesh('findnode z==-5.0');
[Y,I]=sort(FEnode(NodeID,6),'ascend');  NodeID=NodeID(I);
femesh('ObjectBeamLine',NodeID)
femesh('addsel')
%---------------------------------------------------Numeration optimization 
[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));% total nodes
%Title writing
fprintf(fid_geom,'********* TUTORIAL 4 - PARAX 2D **********\n');
%---------------------------------------------------Truncate to 6th decimal
for inod = 1:numnp
    FEnode(inod,5:7)= round(FEnode(inod,5:7)*1e6)*1e-6;
end
%-----------------------------------------------------------Node selections
%Equivalent nodes
NodeID=femesh('findnode y==0.'); NODE_eq=[];
for inode=1:length(NodeID)
    inode_eq=femesh(['findnode z==',num2str(FEnode(NodeID(inode),7)),...
        ' & y==1.0']);
    NODE_eq=[NODE_eq ; NodeID(inode) inode_eq];
end
NODE_eq=[NODE_eq ; NODE_eq(:,2) NODE_eq(:,1)];
%-------------------------------------------------------------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)); %mesh limits
    % id = [x y z rx ry rz p 0 neq]
    % standard node definition
    id=[1 0 0 1 1 1 1 0 0];%2D plane-strain mechanic
    %---------------------------------------------------Boundary conditions
    %Equivalent nodes
    if ismember(inode,NODE_eq(:,1))==1; 
        id(9)=NODE_eq(find(NODE_eq(:,1)==inode),2); end%
    %Lateral limits
    if (abs(FEnode(inode,6)-ymin)<1e-5 || abs(FEnode(inode,6)-ymax)<1e-5) ;
            id(2)=1; end; %fix u_y
    %Write node
    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

%----------------------------------------------------------Elements writing
%NPAR(1)=2: Solid 2D
ITELT=1; IEL=4; IPS=1; EPAI=1.0; MTYP=1; KG=0;%standard parameters
N_group=[20];%number of elements by group
Group_num=[1];%SDT group number
%Quadb elements
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
%NPAR(1)=8: Parax 2D
N_group=1;
group=2;
%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
fprintf(fid_geom,'12345678901234567890123456789012345678901234567890123456789012345678901234567890\n');
%close *.geom file
fclose('all')

En la malla propuesta, hemos dividido la columna de $5$[m]en 20 elementos de $0.25$[m]. En general se asume que se debe contar con $8$ ó $10$ elementos por longitud de onda, de forma que nuestro modelo podrá ''ver'' correctamente ondas de mínimo

$$\lambda \approx 10 \times 0.25 = 2.5 \mbox{ [m]}$$

Como el problema es 1D (y $\nu=0$), la velocidad de propagación de onda $c_p$ queda dada

$$c_p = \sqrt{ \frac{E}{\rho} } = \sqrt{ \frac{10\times10^6}{1600} } \approx 79 \mbox{ [m/s]}$$

Por lo tanto, nuestra malla podrá representar frecuencias de hasta

$$f \leq \frac{79}{2.5} \approx 30 \mbox{ [Hz]}$$

Si se deseara ver ondas de longitud más corta (o frecuencia más alta), se deberá refinar la malla. La malla de elementos finitos generada se ilustra en la Fig.2. La malla cuenta de 20 elementos sólidos 2D, 1 elemento unidimensional (tipo viga) y 42 nodos. El elemento unidimensional permite automatizar la generación de los elementos paraxiales.

Figura 2: Malla empleada

El archivo de comandos para efectuar el cálculo con GEFDyn es el siguiente:

*** TUTORIAL 4: PARAXIAL 2D SOLO ABSORBENTE
   42 1001111 2    1    0    2
    1    2
  150     0.001        0.      1.00
    1    1    0
    1    0    1  -10   11 
   10                                                                               * IOPE
    0                                                                               
    0
         1         1                                                                * IPPGEF/IFORMA
        0.        0.      -10.     1000.        0.
    2    1    1                                                                     * NCHARG/NCOURB/NPTM
         1         3         1        0.         0         0                        * LL/LTYP/NPTC/DATACT/LECHRG/ICHPRI
        1.                                                                          * VALCST
         1         2         1      -.01
         2         2         1      -.01 
    0                                                                               * ICON/IDESEQ
    0                                                                               * JCOUCH
   2  20   1   0   1   2   4   0   0   2   0   0   1   0   1   1           0        * NPAR(1)=2, NPAR(8)=0, NPAR(15)=1 
    1     1600.                                                                     * N/DEN/POROS0
    10.E06        0.        0.    20.e10        0.        1.                        * E/XNU/PHI/C/PSI/AK0
        0.        0.         2                                                      * FACNAP/ZNAPPE/IMODEL
        0                                                                           * NCOUCH
         1                                                                          * IPPGEF
   8   1   1   0   1   2   0   0   0   2   0   0   0   0   1   1               1    * NPAR(1)=8, NPAR(14)=0 
    1     1600.
    10.0E6       0.0
         1         0         0
123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890

Sólo describiremos en detalle algunos aspectos particulares del problema. El significado de los otros campos puede ser interpretado de acuerdo a la información del manual.

El primer aspecto relevante corresponde al código de integración en tiempo. En este caso se ha empleado la opción IOPE=10 ya que se requiere de un esquema explícito/implícito. Los elementos de volumen 2D son integrados implícitamente (podrían ser explícitos), sin embargo los elementos paraxiales siempre son integrados explícitamente. Ya que en el modelo hay ambos tipos de elementos, se debe emplear el código 10. Es importante notar que son justamente los elementos paraxiales por el hecho de ser integrado explícitamente los que controlan el paso de integración. En este caso se ha empleado $\Delta t =0.001$ que es el valor usual para casos sísmicos. Si el mallado fuera más fino probablemente habría que emplear un paso también más fino. Se deja al lector chequear que para un paso un poco más grande el cálculo diverge.

Para imponer el esfuerzo $\sigma_0$ en superficie, impondremos una fuerza nodal de 0.01 [N] (hacia abajo), como el ancho de la columna es unitario, esta fuerza impone una compresión de $\sigma_0=0.01 \times 2 / 1 = 0.02$ [Pa]. Como el valor es constante, se empleó una curva de tipo LTYP=3.

Se imponen condiciones nodales nulas (ICON=0) y como no se consideran fuerzas de gravedad no se requiere inicializar los esfuerzos internos (IDESEQ=0).

El suelo se modelo mediante 20 elementos sólidos 2D (NPAR(1)=2), no se considera la gravedad (NPAR(8)=0) y asumiremos comportamiento de tipo Mohr-Coulomb (NPAR(15)=1). Como se asume material elástico, usaremos un valor muy grande para la cohesión (C=20.e10).

La imposición del borde absorbente interior se efectuará mediante elementos paraxiales 2D (NPAR(1)=8). En este ejemplo emplearemos estos elementos sólo para absorber la onda inducida por la carga en superficie, de forma que NPAR(14)=0. Como el elemento debe reemplazar un medio infinito, sus propiedades materiales deben coincidir con las del suelo ($\rho$, $E$ y $\nu$). En el caso de esto elementos, es fundamental su orientación. De acuerdo a lo indicado en la descripción de los elementos, la numeración de los elementos debe inducir una normal hacia el interior del medio (numeración de izquierda a derecha en este caso).