How to input variable mechanical properties for the structural mechanics solver?

Hi every one,

I've been trying to get my head around Matlab's FEM structural mechanics solver (PDE toolbox), and there is one issue I'm still struggling with:

Is it possible to input mechanical properties (mainly the Youngs modulus) which would vary with temperature / location? The structuralProperties function will only take scalar values, but I'd like to feed it a function handle to account for space-dependent properties, i.e something like E = @(location,state) E0 * location.x.*location.y .

I am running the structural model on a structure with a (known) non-uniform temperature profile. Mechanical properties being temperature dependent, they will therefore have different values in different parts of my structure.

Thank you for your help,

[Patrick Frey]

Is there really no solution???

[Tejas Ladhe]

Did you get any solution? I am stuck in the same problem in which I want to define structural properties as function of temperature (state.u). Plz help.

[Thibaut Desguers]

No, I've ended up coding it myself...

[alexander fernandez]

i have this same problem,are you able to share you code please?

Answer 1

[Roman Kandinskii]

Hello Thibaut,

It's been a while since your original post, but I have figured out how to make the PDE toolbox accept the Young's modulus as a function of coordinates.

It was a nasty ride, so if there's anyone interested, buckle up.

As noted by Thibaut, the 'structural' model object accepts the material properties only as constants. Who knows what difficulty it was to make them definable as functions (in the thermal objects you can set a variable conductivity easily). So the key is to use the general PDE object as described in, for example, here. The general PDEmodel object is more flexible than dedicated structural or thermal objects - you can define your PDE(s) the way you want and specify all the coefficients the way you want. The drawback, one of them at least, is that you have to interpret all the results yourself. Recalculate displacements into stresses, etc.

Let's look at an example. I want to stretch a square plate with a cylindrical inclusion of higher stiffness. You can see it as a composite unite cell (a fibre in matrix).

model = createpde(3);                         % 3 variables - ux, uy and uz. This example is in 3D
g = importGeometry(model,'square_flat.stl');  % using an imported geometry is easier in 3D. The .stl file is attached
extrude(g,10);                              % extrudes the square into parallelepiped
pdegplot(model,'EdgeLabel','on','VertexLabels','on','FaceLabels','on');

Generate the mesh

Hmax = 3.5;
msh = generateMesh(model,'GeometricOrder','quadratic','Hmax', Hmax);
figure;  pdemesh(model);
title('Mesh with Quadratic Triangular Elements');
dofss = length(msh.Nodes)              % Number of degrees of dreedom

Boundary conditions

applyBoundaryCondition(model,'dirichlet','Face',3,'u',[0,0,0]);    % fix one side
applyBoundaryCondition(model,'dirichlet','Face',5,'u',[110,0,0]);  % apply a displacement to the opposite side

Here is a trick to define constant isotropic material for the general PDE object. I found it in the depths of the documentation. It uses an undocumented function elasticityC3D, inside which is this:

function cmat = elasticityC3D(E,nu)
% ELASTICITYC3D Calculate c-matrix for isotropic elastic solid
%   This convenience function is included in the toolbox for the use in
%   several examples. This undocumented function may be removed in a future
%   release.
%
%   This function uses the linearized small-displacement assumption for an
%   isotropic material. Hence, there are only two independent elasticity
%   constants, Young's modulus (E) and Poisson's ratio (nu), which are used
%   to construct c-matrix in PDE Toolbox format.

%   Copyright 2014-2016 The MathWorks, Inc.
G = E/(2*(1+nu));
c1 = E*(1-nu)/((1+nu)*(1-2*nu));
c12 = c1*nu/(1-nu);
C11 = [c1 0 G 0 0 G];
C12 = [0 G 0 c12 0 0 0 0 0];
C13 = [0 0 G 0 0 0 c12 0 0];
C22 = [G 0 c1 0 0 G];
C23 = [0 0 0 0 0 G 0 c12 0];
C33 = [G 0 G 0 0 c1];
cmat  = [C11 C12 C22 C13 C23 C33]';
end

SO, this function takes the Young's modulus and Poisson's ratio, and assembles the c-coefficient matrix for the PDE system in the required vectorial format. It is indeed the stiffness matrix of the elasticity equation, which is described there https://nl.mathworks.com/help/pde/ug/3-d-linear-elasticity-equations-in-toolbox-form.html

The specificity of the vectorial form of the c-coefficient matrix is described here: https://nl.mathworks.com/help/pde/ug/c-coefficient-for-systems-for-specifycoefficients.html#bu5xabm-21 . It is pretty hectic, but it gives the infromation on how to set this matrix as spatially dependent. I doubt if I could wrap my head around it, if not the elasticityC3D function.

If we opt for the constant modulus, we define the PDE coefficients as follows:

E = 200e9;
nu = 0.3;
specifyCoefficients(model,'m',0,...
                          'd',0,...
                          'c',elasticityC3D(E,nu),...
                          'a',0,...
                          'f',[0;0;0]);

Solve the system:

rslt = solvepde(model);
u = rslt.NodalSolution;
pdeplot3D(model,'ColorMapData',u(:,1)); title('ux')
pdeplot3D(model,'ColorMapData',u(:,2)); title('uy')
pdeplot3D(model,'ColorMapData',u(:,3)); title('uz')

And we see the displecement fields as expected

Now, instead of using the elasticityC3D, let's define our own function, where the non-zero coefficient matrix entries are going to depend on coordinates. I am using a sort of logistic function to draw a circular field with certain slope. But it does not necessarily have to be like this.

function cmatrix  = ccoeffunction(location,state)

% cylinder logistic function
E_m = 100E2;    % modulus of the matrix
E_f = 10E9;     % modulus of the fibre
nu = 0.3;       % Poisson's ratio
kk = 3000;      % steepness of the logistic curve
r = 35;         % radius of the fibre
x0 = 50; y0 = 50;  % centre coordinates

% The logistic function. I am using dot operators (.* ./ .^) everywhere in this function
% to ensure element-wise operations
E = E_m + (E_f -E_m)./(  1 + exp(  kk*(  sqrt((location.x-x0).^2+(location.y-y0).^2) - 1*r))   ); 

% the vectorial form of the matrix is now a matrix (yeah, I know right), of size n-by-m, 
% where n is the number of entries in our particular form of the c-coefficient matrix,
% and m is the number of degrees-of-freedom (number of mech nodes)
cmatrix  = zeros(45,numel(location.x));

G = E./(2*(1+nu));                  % again, use dots to avoid matrix operations
c1 = E.*(1-nu)./((1+nu)*(1-2*nu));
c12 = c1.*nu/(1-nu);

cmatrix (1,:) = c1;
cmatrix (3,:) = G;
cmatrix (6,:) = G;
cmatrix (8,:) = G;
cmatrix (10,:) = c12;
cmatrix (16,:) = G;
cmatrix (18,:) = c1;
cmatrix (21,:) = G;
cmatrix (24,:) = G;
cmatrix (28,:) = c12;
cmatrix (36,:) = G;
cmatrix (38,:) = c12;
cmatrix (40,:) = G;
cmatrix (42,:) = G;
cmatrix (45,:) = c1;
end

How exctly the "vectorial" from relates to the c-coefficient matrix, you can find here https://nl.mathworks.com/help/pde/ug/c-coefficient-for-systems-for-specifycoefficients.html#bu5xabm-10

See the section 3N(3N+1)/2-Element Column Vector c, 3-D Systems.

Let's now re-evaluate the coefficients and solve the problem.

specifyCoefficients(model,'m',0,...
                          'd',0,...
                          'c',@ccoeffunction,...
                          'a',0,...
                          'f',[0;0;0]);
rslt = solvepde(model);
u = rslt.NodalSolution;
pdeplot3D(model,'ColorMapData',u(:,1)); title('ux')
pdeplot3D(model,'ColorMapData',u(:,2)); title('uy')
pdeplot3D(model,'ColorMapData',u(:,3)); title('uz')

And now we see the fictitious fibre in the displacement fields.

During the solve, there was a warning of poor conditioning. So it's probably better to define different geometric domains for such problem. But this suffices for the example.

I hope it helps someone,

Roman

[Roman Kandinskii]

A little comment: numel(location.x) is not exactly the number of degrees of freedom, so be sure you use it as the second size of your c-coefficient matrix.

[Spencer Dansereau]

Awesome modification! I'd like to use your method here, but adapt it to what I'm doing. I'm trying to have not just a fiber in the middle, but adapt the whole volume to a 2D map of high and low E values. Basically I have a 100x100 matrix of 1's and 0's, that I'd like to map onto the 100x100x10 block, with each element within the mesh being mapped to the corresponding E value dictated by my image map. However I'm running into a problem where the E value for each time the function is run in the solver is not necessaru a 1x1 value, but rather is occationally a 1x2500 value. I'm confused on why this is first off, but also it causes a problem where I can't directly map the E values from the image map onto the coefficents.Getting the image into the function is just through a global variable I set as the image, which also works with the function you wrote with the x0, y0 and r values if you want to change that.