2 Generación de la malla

Para construir una malla de elementos finitos para el problema emplearemos las herramientas de SDT. Por supuesto, esta etapa puede ser efectuada con cualquier programa de generación de mallas, generando el archivo de geometría en el formato de lectura de GEFDyn.

Comenzaremos definiendo el tamaño del dominio del problema. En primer lugar, dada la simetría, podemos emplear un mallado correspondiente sólo a la mitad del dominio. En este caso emplearemos la mitad derecha, luego la fundación quedará representada por un ancho de 1.0 [m]. Como primera aproximación fijaremos la extensión del malla como 10 veces en el sentido horizontal. Las condiciones de borde quedarán definidas por un desplazamiento horizontal nulo en los bordes laterales, mientras que el borde inferior queda completamente limitado en el sentido vertical.

Como primera aproximación emplearemos una malla relativamente gruesa y homogénea, constituida por elementos cuadrángulos lineales de $0.5 \times 0.5$[m] de costado. La malla empleada se muestra en la Figura 2.

Figura 2: Malla empleada

Para construir la malla con SDT se empleó el código siguiente:

filename='tut2_footing';%filename
FEnode=[];  FEelt=[]; %reseting default variables
%---------------------------------------------------------Mesh construction
FEel0=[]; 
FEnode=[1 0 0 0 0.  0 .0
        2 0 0 0 0. 10 .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 -4')
femesh('divide 8 20')
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 2 - RIGID FOOTING **********\n');
%---------------------------------------------------Truncate to 6th decimal
for inod = 1:numnp
    FEnode(inod,5:7)= round(FEnode(inod,5:7)*1e6)*1e-6;
end
%-----------------------------------------------------------Node selections
Node_foot=femesh('findnode z==0. & y<=1.0'); %footing nodes
%-------------------------------------------------------------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
    %Footing nodes
    if ismember(inode,Node_foot)==1; id(3)=-1; end%imposed vertical DOF 
    %Lateral limits
    if (abs(FEnode(inode,6)-ymin)<1e-5 || abs(FEnode(inode,6)-ymax)<1e-5) ;
            id(2)=1; end; %fix u_y
    %Bottom limit
    if abs(FEnode(inode,7)-zmin)<1e-5 ;
            id(3)=1; end; %fix u_z
    %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=[160];%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
fprintf(fid_geom,'12345678901234567890123456789012345678901234567890123456789012345678901234567890\n');
%close *.geom file
fclose('all')

La primera parte del código (Mesh construction) define la construcción de la malla empleando funciones de SDT. Consultar las opciones de femesh para mayores detalles.

La sección siguiente (Numeration optimization) permite optimizar la enumeración de los nodos de forma de minimizar el ancho de banda de la matriz de rigidez tangente. Esta opción es importante cuando se emplea la versión secuencial de GEFDyn, que emplea el sistema de almacenaje skyline. En la versión paralela del código emplea matrices sparse de forma que la enumeración de los nodos es irrelevante. En dicho caso es más importante la generación de las particiones de la malla.

Luego, se escribe el archivo *.geom en el formato adecuado para GEFDyn. Las líneas siguientes crean el archivo, lo abren y escriben el encabezado. Para evitar problemas de truncado de dígitos al crear el archivo de geometría, se limita el número de posiciones decimales a $6$ en la línea siguiente.

Para definir las condiciones de borde, es frecuente seleccionar grupos de nodos cuyas condiciones son idénticas. En este caso se seleccionaron los nodos que tendrán un desplazamiento vertical impuesto vertical.

A continuación se procede a la escritura de los nodos mediante un loop sobre la definición del mallado. En este archivo de ejemplo, se ha decidido emplear una definición estándar de las condiciones iniciales:

    % standard node definition
    id=[1 0 0 1 1 1 1 0 0];%2D plane-strain mechanic

De acuerdo a los tipos de problemas de GEFDyn, en un problema mecánico 2D en deformaciones planas sólo se liberan los desplazamientos horizontales $u_y$ y los verticales $u_z$ (2da y 3ra posición de id). En efecto el vector id posee 9 componentes: las primeras $7$ definen el código de los grados de libertad, mientras que las dos últimas se emplean para definir condiciones de equivalencia (no empleadas en este ejemplo).

Previamente a la definición de id, se calcularon las coordenadas máximas y mínimas segun $y$ y $z$. Estos valores permiten imponer las condiciones de borde. Por ejemplo, se bloquea el borde inferior id(3)=1 cuando la coordenada $z$ del nodo es igual a la mínima zmin. La imposición de las otras condiciones de borde se deduce directamente del código Matlab.

Finalmente se procede a la escritura de los tipos de elementos. En general se suele emplear loop para grupos elementos topológicamente idénticos. En el caso de grupos de tipo distintos, se pueden ir definiendo uno a continuación del otro.

En este caso existen un sólo tipo de grupo (volumen 2D) compuestos por $160$ elementos. El orden de escritura en el archivo *.geom define el número del grupo en GEFDyn, pero no tiene porqué ser el mismo definido en la construcción de la malla SDT. La variable Group_num permite indicar dicha correspondencia. Se define también los valores correspondientes a una serie de parámetros de este tipo de elementos (volumen 2D) para GEFDyn. Luego, los elementos se escriben uno a uno mediante un loop con el formato adecuado.

Finalmente el archivo se cierra y hemos completado la creación de la malla para GEFDyn.

Observación: durante la construcción de la malla existen varios comandos que nos permiten controlar/verificar si está consiguiendo el resultado deseado:

• femesh plotelt traza la malla actualmente asociada a FEelt y FEnode.
• femesh plotel0 traza la malla actualmente asociada a FEel0 y FEnode (selección de trabajo)
• femesh info despliega en pantalla el tipo y cantidad de elementos contenidos en el modelo (FEelt), en la selección de trabajo (FEel0) y el número de nodos definidos en FEnode.