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.