Refining mesh size (by having very small mesh size) close to the semi-circular flaw tip - Maximum Principal Stresses and Maximum Shear Stresses

Good morning,

I hope this e-mail finds you well. Your input will be of incredible importance in my thesis. I am using the PDE (Partial Differential Equation) ToolBox of Matlab to obtain some stress distributions. Basically regarding the geometry, I am modeling a flaw (opening with the geometry of interest visualized upon running the below code). The major problem with the PDE ToolBox is the meshing. I am basically interested in obtaining the values of stress close to the flaw tip (tip of the opening), Therefore, I want my mesh size close to that region to be very very very small. So far, on my laptop, I am able to use a maximum mesh size (Hmax) equal to 0.135mm (smaller than that, it does not work) but i definitely need to get much smaller than this to have precise results close to the flaw tip. Because, with the actual mesh size that I have, I am not obtaining precise results. Your input here will be of extreme importance. You can find enclosed my code and the biggest challenge/problem as mentioned before is the mesh size. I am very interested in having an extremely small mesh size close to the flaw tip.

CODE:

___________________________

model=createpde('structural','static-planestress');

%flaw width b taken to be 0.032";
b=0.032*25.4;

%Radius of the flaw tip equivalent to half of the flaw width b;
Radius=b/2;

%L: Length of the flaw;
L=0.5*25.4;

%La: Length of the rectangle;
La=L-b;

xsemiL = 1*25.4;
ysemiL = 2*25.4;

% semiL= 50;

alpha = 10; 

R = [3 4 -xsemiL xsemiL xsemiL -xsemiL -ysemiL -ysemiL ysemiL ysemiL]';

R1 = [3 4 (La/2)*cosd(alpha)-Radius*cosd(90-alpha) (La/2)*cosd(alpha)+Radius*cosd(90-alpha) -(La/2)*cosd(alpha)+Radius*cosd(90-alpha) (-La/2)*cosd(alpha)-Radius*cosd(90-alpha) ...
    (La/2)*sind(alpha)+Radius*sind(90-alpha)  (La/2)*sind(alpha)-Radius*sind(90-alpha) -(La/2)*sind(alpha)-Radius*sind(90-alpha) -(La/2)*sind(alpha)+Radius*sind(90-alpha)]';

C1 = [1 (La/2)*cosd(alpha) (La/2)*sind(alpha) Radius 0 0 0 0 0 0]';

C2 = [1 (-La/2)*cosd(alpha) (-La/2)*sind(alpha) Radius 0 0 0 0 0 0]';

gdm = [R C1 C2 R1];
ns = char('R','C1','C2','R1');
sf='R-C1-C2-R1';
g = decsg(gdm,sf,ns');

geometryFromEdges(model,g);

figure
pdegplot(model,'EdgeLabel','on','VertexLabel','on', 'Facelabel','on');
axis equal 
title 'Geometry with Edge, Face and Vertex Labels';

structuralProperties(model,'Face',1,'YoungsModulus',30E3,'PoissonsRatio',0.25,'MassDensity',2.63E-6);

structuralBC(model,'Edge',5,'YDisplacement',0);
structuralBC(model,'Vertex',4,'XDisplacement',0);

structuralBoundaryLoad(model,'Edge',2,'SurfaceTraction',[0;-5]);

generateMesh(model, 'Hmax', Radius/5)

figure
pdemesh(model);

___________________________

Please, keep me informed if you need any additional information. I can't wait to receive from you.

Best regards,

Mohamad Zaarour

Answer 1

[Alan Weiss]

I really don't know, but it is possible that the legacy function adaptmesh would enable you to perform your calculation to the requisite accuracy. You cannot use a pde model, and you'd have to write your own equations, so the setup is awkward, but the result might be worthwhile.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

Answer 2

[Konstantin Kovalev]

Hello Mohamad,

Thank you for interest in PDE Toolbox. I may have a suggestion for you that could work. As I understand, you are modeling a crack, want to capture extreme stresses at tip, and thus need high resolution mesh there that you can't get with the toolbox.

I am not sure how important is that the tip shape is circular but if you have some flexibility there, you might consider using elliptical shape instead (see the image). In this case the ellipse's eccentricity can be used to control the max curvature and thus the mesh size because curvature is the main mesh size driver (make sure to specify small enough min size).

Please note that this is a suggestion, I did not write code to try it.

I hope this helps. Let me know if it doesn't work for you.

Thanks,

Konstantin