I am trying to put together a basic script to generate a 2D FEA stress analysis for a thin plate with various configurations of holes(by different diameters, locations, and # of holes) so I can narrow down the configurations before going to a full 3D FEM in Abaqus or something that allows for more advanced analyses. I am following an example from the MATLAB PDE documentation, https://www.mathworks.com/help/pde/ug/stress-concentration-in-plate-with-circular-hole.html, and my geometry and mesh seems to generate as I expect but I am running into an issue with my BCs and applying loads.
In the example they use a surface traction with arguments [100;0] but what exactly these values represent is not specified in the example or the documentation, https://www.mathworks.com/help/pde/ug/edgeload.html, and for some reason the surface traction here seems to include a normal distributed force even though a surface traction usually exclusively refers to tangential forces. I am not sure what units is intended by the [100;0] arguments.
For the example, they apply a distributed force normal to the right edge and constrain the x-displacement of the right edge and the y-displacement of the bottom left vertex. I believe my force loading is fairly similar (uniaxial tensile stress uniformly distributed along the horizontal edges, and therefore I tried to use similar BCs as the example however my FEA stress results seem to produce something very different. There is a noticeable contour color change when I zoom in on the bottom right corner where I applied my vertex BC though, which I'm not sure the cause of. Should I use a pressure instead for applying my load or is there a better way to go about this issue? Below are some figures relevant to my problem formulation and code results.
Problem setup:
Code:
% define plate geometry
% plate dimensions
halfWidth = 125;
halfLength = 125;
numberOfHoles = 5;
% hole radii and locations
[radius1, x_circle1, y_circle1] = deal(10/2, 0, 0);
[radius2, x_circle2, y_circle2] = deal(20/2, 50, 50);
[radius3, x_circle3, y_circle3] = deal(20/2, -50, 50);
[radius4, x_circle4, y_circle4] = deal(20/2, -50, -50);
[radius5, x_circle5, y_circle5] = deal(20/2, 50, -50);
% define the geometry description matrix (GDM)
% rectangle
R1 = [3 4 -halfLength halfLength ...
halfLength -halfLength ...
-halfWidth -halfWidth halfWidth halfWidth]';
% circle
C1 = [1 x_circle1 y_circle1 radius1 0 0 0 0 0 0]';
C2 = [1 x_circle2 y_circle2 radius2 0 0 0 0 0 0]';
C3 = [1 x_circle3 y_circle3 radius3 0 0 0 0 0 0]';
C4 = [1 x_circle4 y_circle4 radius4 0 0 0 0 0 0]';
C5 = [1 x_circle5 y_circle5 radius5 0 0 0 0 0 0]';
% define the combined GDM, name-space matrix, and set
% formula to construct decomposed geometry using decsg
if numberOfHoles == 3
gdm = [R1 C1 C2 C3];
ns = char('R1','C1','C2','C3');
g = decsg(gdm,'R1 - C1 - C2 - C3',ns');
elseif numberOfHoles == 4
gdm = [R1 C1 C2 C3 C4];
ns = char('R1','C1','C2','C3','C4');
g = decsg(gdm,'R1 - C1 - C2 - C3 - C4',ns');
elseif numberOfHoles == 5
gdm = [R1 C1 C2 C3 C4 C5];
ns = char('R1','C1','C2','C3','C4','C5');
g = decsg(gdm,'R1 - C1 - C2 - C3 - C4 - C5',ns');
end
figure
pdegplot(g,EdgeLabel="on");
axis([-1.2*halfLength 1.2*halfLength -1.2*halfWidth 1.2*halfWidth])
title("Geometry with Edge Labels")
figure
pdegplot(g,VertexLabels="on");
axis([-1.2*halfLength 1.2*halfLength -1.2*halfWidth 1.2*halfWidth])
title("Geometry with Vertex Labels")
model = femodel(AnalysisType="structuralStatic", Geometry=g);
% Properties and loading
model.MaterialProperties = materialProperties(YoungsModulus=200E3, PoissonsRatio=0.25);
model.EdgeLoad(2) = edgeLoad(SurfaceTraction=[-100;0]); % negative for tension
% BCs
model.EdgeBC(4) = edgeBC(YDisplacement=0);
model.VertexBC(1) = vertexBC(XDisplacement=0);
model_mesh = generateMesh(model,Hmax=5/6);
figure
pdemesh(model_mesh)
title('Model Mesh')
R = solve(model_mesh);
figure
pdeplot(R.Mesh,XYData=R.Stress.sxx, ColorMap="jet")
axis equal
title("Normal Stress Along y-Direction")
Bottom right corner contour results(where the vertex BC of is applied):
