Dividing output of structuralmodel solution based on their cell.

I'm trying to do 3D modelling of stacked cylinder under some static pressure using structuralModel. I currently have 2 cylinder, where on is stacked on top of another and applying some static pressure on the uppermost surface. I solved the model and trying to get the stress and strain from each part of the mesh. I don't need exact coordinates of the values but need to classify which cell (upper or lower part of cylinder) that they come from. Below are my codes. Don't worry about the HGforce function as its just function that I use to create some static pressure profile.

in = 25.4e-3

structuralModel = createpde("structural","static-solid");

gm = multicylinder(25.4e-3,[0.25*in,266e-9], "Zoffset",[0,0.25*in]);

generateMesh(structuralModel, 'Hmax', 2.0e-3);

E1 = 6.8e10;

sigma1 = 0.18;

E2 = 6.9e10;

sigma2 = 0.19;

structuralProperties(structuralModel, "Cell", 1, ...

'YoungsModulus', E1, ...

'PoissonsRatio', sigma1);

structuralProperties(structuralModel, "Cell", 2, ...

'YoungsModulus', E2, ...

'PoissonsRatio', sigma2);

% Set back face of mirror fixed boundary condition

structuralBC(structuralModel, 'Face', 1, 'Constraint', 'fixed'); % Fix back face of mirror

structFunForce = @(location,state) HGforce(location, state, 0, 0, F0, w0);

structuralBoundaryLoad(structuralModel, ...

'Face', 4, ...

'SurfaceTraction', structFunForce);

R = solve(structuralModel);

Answer 1

[Paras Gupta]

Hello,

I understand that you want to get the stress and strain for each mesh part of structuralmodel into the upper or the lower cell and classify the part into the cell it belongs to. I assume parts of mesh mean the Elements in a FEMesh object. Please find below the code to achieve the same.

in = 25.4e-3;
F0 = 100; % Assuming a value for F0
w0 = 1; % Assuming a value for w0

structuralModel = createpde('structural','static-solid');

cylinderGeom = multicylinder(25.4e-3,[0.25*in,0.25*in], "Zoffset",[0,0.25*25.4e-3]);

% Assign the geometry to the structural model
structuralModel.Geometry = cylinderGeom;

% Set the mesh size
generateMesh(structuralModel, 'Hmax', 2.0e-3);

E1 = 6.8e10; 
sigma1 = 0.18; 
E2 = 6.9e10; 
sigma2 = 0.19; 
structuralProperties(structuralModel, "Cell", 1, ...
    'YoungsModulus', E1, ...
    'PoissonsRatio', sigma1);
structuralProperties(structuralModel, "Cell", 2, ...
    'YoungsModulus', E2, ...
    'PoissonsRatio', sigma2);
structuralBC(structuralModel, 'Face', 1, 'Constraint', 'fixed'); % Fix back face of mirror
structFunForce = @(location,state) HGforce(location, state, 0, 0, F0, w0);
structuralBoundaryLoad(structuralModel, ...
        'Face', 4, ...
        'SurfaceTraction', structFunForce);

structuralResults = solve(structuralModel);

% Cell arrays to store szz component of stress and index of every mesh
% element based on whether they belong to upper cylindrical cell or lower
upperElements = {};
lowerElements = {};

% % Loop over each element and classify based on z-coordinate
numElements = size(structuralModel.Mesh.Elements, 2);

for i = 1:numElements
    % For each element, find the corresponding node indices
    elementNodes = structuralModel.Mesh.Elements(:, i);
    % Find the coordinates of each node in the element
    elementCoordinates = structuralModel.Mesh.Nodes(:, elementNodes);
    % Find the szz component of stress for the element under consideration
    stress_element_z = mean(structuralResults.Stress.szz(elementNodes));
    % Find the centroid of the current element by taking mean
    centroid = mean(elementCoordinates, 2);
    % If the z-coordinate of the centroid is above 0.25*in, we classify the
    % element to come from upper cell, and vice versa
    zCoord = centroid(3);
    if zCoord >= 0.25*in
        upperElements{end+1} = {i, stress_element_z};
    else
        lowerElements{end+1} = {i, stress_element_z};    
    end
end

% Dummy HGForce
function force = HGforce(location, state, ~, ~, F0, w0)
    force = [0; 0; F0 * cos(w0 * state.time)];
end

The above code finds the szz component of stress (FEStruct object) for each mesh element. We can similarly find the other components of stress and repeat the process for strain (FeStruct Object).

You can refer to the following documentations for more information on the code above:

• https://www.mathworks.com/help/pde/ug/pde.femesh-properties.html
• https://www.mathworks.com/help/pde/ug/pde.staticstructuralresults.html

Hope this helps.