How can I apply boundary conditions to a random point in my Geometry for my FEM in PDE Toolbox (Not just define BCs at a face, edge, vertex)?

Hi, I am trying to simulate a simple 3 point bending test using the PDE Toolbox. I'm struggling with getting the correct boundary conditions. My test piece is a rectangular prism with a load in the -z direction on the top face, and since it is 3 point bending it should be supported at two points which are not allowed to move in the -z. These points should be 30 mm in from each of the side faces of the prism in the axial direction. However I haven't been able to figure out how to apply boundary conditions to arbitrary points of my geometry. I only see the options to apply them to faces, edges, and vertices, which doesn't help with 3 point bending. My code is below, thanks for help in advance.

clear all

close all

tic % times total time of script

% Define Material Properties

E = 70E9; % Young's Modulus in Pa

nu = 0.34; % Poisson's Ratio

load = 2000;

distributed_load = 2*10^7; % 2kN load distributed over the face

depth_mm = 1; % Depth in mm

width_mm = 5; % Width in mm

% Script for FEM using PDE Toolbox

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

% Import geometry here from an STL file - make it

importGeometry(model,"PerpendicularCase.stl")

% Plot the Geometry and Label it

pdegplot(model,'FaceLabels','on', 'VertexLabels', 'off', 'FaceAlpha',0.5); % Alpha is transparency value, face labels labels things

% Discretize geometry into finite elements

msh = generateMesh(model,'Hmax',1); % Hmax is max edge length of elements. One is a good amount in this case

%pdeplot3D(model)

%%

% Create Physics of Model

structuralProperties(model, 'YoungsModulus', E, 'PoissonsRatio', nu);

%Define Loads

%Ignoring Clamp Length for Now

%Top face is F10

%Bottom Face is F3 - side faces are F2 and F4, want to apply BC 30 mm in from them

structuralBoundaryLoad(model,'Face',10,'SurfaceTraction',[0,0,-distributed_load])

result = solve(model);

% displace = result.Displacement;

% z_displace = displace.uz;

%

% max_abs_z_displacement = max(abs(z_displace));

%

% fprintf('Maximum z displacement: %f\n', max_abs_z_displacement);

% Plot solution

pdeplot3D(model, 'ColorMapData', result.VonMisesStress, 'Deformation', result.Displacement, 'DeformationScaleFactor', 1);

title(sprintf("Von Mises Stresses for depth = %dmm, width = %dmm", depth_mm, width_mm))

average = mean(result.VonMisesStress);

max_stress = max(result.VonMisesStress);

fprintf('The average Von Mises Stress is: %.2f MPa\n', average/(10^6));

fprintf('The maximum Von Mises Stress is: %.2f MPa\n', max_stress/(10^6));

maxdisplacement = max(abs(result.Displacement.z));

fprintf('The maximum Z displacement is: %f mm\n', -maxdisplacement);

%Sort out extreme stresses

threshold = 10^10;

filtered_stress = result.VonMisesStress(result.VonMisesStress<=threshold);

adjusted_average = mean(filtered_stress);

adjusted_max = max(filtered_stress);

toc

Answer 1

[Abhishek Chakram]

Hi Andrew Olson,

Applying boundary conditions to specific points (rather than to entire faces, edges, or vertices) in the geometry for FEM analysis in PDE Toolbox can be challenging because the toolbox is designed to work with boundary conditions applied to geometric entities. However, for simulations such as a three-point bending test, where you want to apply loads or constraints at specific points, you can use a workaround by creating a very small region around those points and applying the boundary conditions to that region.

Since your geometry is imported from an STL file, consider modifying the file to include small cylindrical or spherical regions at the points where the supports are located. This might require using 3D CAD software before importing the geometry into MATLAB. These regions can then be treated as separate entities (faces) in the FEM model.

After modifying your geometry and importing it into MATLAB, you can apply boundary conditions to the faces corresponding to the small regions you added.

Here is example code for the same:

clear all
close all
tic % times total time of script

% Define Material Properties
E = 70E9; % Young's Modulus in Pa
nu = 0.34; % Poisson's Ratio
distributed_load = 2*10^7; % 2kN load distributed over the top face
depth_mm = 1; % Depth in mm
width_mm = 5; % Width in mm

% Script for FEM using PDE Toolbox
model = createpde('structural','static-solid');

% Import geometry here from an STL file - make it
importGeometry(model,"PerpendicularCase.stl");

% Plot the Geometry and Label it
pdegplot(model,'FaceLabels','on', 'VertexLabels', 'off', 'FaceAlpha',0.5); % Alpha is transparency value, face labels labels things

% Discretize geometry into finite elements
msh = generateMesh(model,'Hmax',1); % Hmax is max edge length of elements. One is a good amount in this case

% Create Physics of Model
structuralProperties(model, 'YoungsModulus', E, 'PoissonsRatio', nu);

% Define Loads
% Assuming top face is labeled as 10 in your geometry
structuralBoundaryLoad(model,'Face',10,'SurfaceTraction',[0,0,-distributed_load]);

% Assuming you have identified the faces or created small regions for supports, apply boundary conditions
% For example, if the small support regions are faces 11 and 12 in your geometry
structuralBC(model,'Face',11,'Constraint','fixed');
structuralBC(model,'Face',12,'Constraint','fixed');

% Solve the model
result = solve(model);

% Plot solution
pdeplot3D(model, 'ColorMapData', result.VonMisesStress, 'Deformation', result.Displacement, 'DeformationScaleFactor', 1);
title(sprintf("Von Mises Stresses for depth = %dmm, width = %dmm", depth_mm, width_mm));

% Calculations and Output
average = mean(result.VonMisesStress);
max_stress = max(result.VonMisesStress);
fprintf('The average Von Mises Stress is: %.2f MPa\n', average/(10^6));
fprintf('The maximum Von Mises Stress is: %.2f MPa\n', max_stress/(10^6));
maxdisplacement = max(abs(result.Displacement.z));
fprintf('The maximum Z displacement is: %f mm\n', -maxdisplacement);

toc

You can refer to the following documentation to know more about functions used:

• structuralBC: https://www.mathworks.com/help/pde/ug/pde.structuralmodel.structuralbc.html

Best Regards,

Abhishek Chakram

[Andrew Olson]

I figured I had to do something like that, thank you for your help!