Definition of the c coefficient for assempde
3 次查看(过去 30 天)
显示 更早的评论
Hello, I am working on an elasticity problem and I would like to solve the mechanical equilibrium equation with the pde toolbox (f=0, a=0). I work in 2D and have periodic boundary conditions. My definition of the inhomogeneous stiffness tensor (c in the pde) depends on a function defined in the whole domain (depending on x,y). The inhomogeneous stiffness tensor can be described through a 4x4 matrix (entries C_1111, C_1112, C1121...depending on the position). I have already tried to define the matrix as indicated in the manual, but also reducing my problem into a 3x3 matrix, my entries have dimensions of (nr. of triangles in the mesh (times) nr. of triangles in the mesh).. Can I use the pde toolbox to solve my problem? Thanks in advance!
Laura
0 个评论
采纳的回答
Bill Greene
2013-12-3
Hi,
Yes, you can solve your problem with PDE Toolbox.
I recommend writing a MATLAB function to define your c-coefficient as shown here for the scalar case:
For your elasticity case, as you note, your c-coefficient can be thought of as a 4x4 matrix. There are several ways to define the terms in this matrix. If you want to define all the entries in the upper triangle, you would need to define 10 terms. As you note, the MATLAB function would need to return a matrix with dimensions 10 x number_of_triangles. Since your c-matrix is a function of x and y, each column in that matrix could have different values. The documentation page I list above shows how to calculate positions of the element centroids, solution values, and gradients that you might need to define your coefficients as a function of position.
Bill
3 个评论
Bill Greene
2013-12-4
Hi,
You can pass a 10 x number_of_triangles sized c-matrix to assempde like this assempde(b,p,e,t,c,...) Each of the terms in this matrix must be a scalar.
It sounds like you are running into difficulties computing the terms in your c-matrix from the results of your eta(x,y) function. If you post the code you use for this step, we may be able to suggest the appropriate code so you get a 10 x number_of_triangles matrix at the end of this step.
Bill
更多回答(2 个)
Bill Greene
2013-12-5
Hi,
Your coefficients are sufficiently complicated that I'm not sure I understand all the issues, but hopefully I can at least point you in a direction where you can resolve the details, yourself.
Looking at your code snippet, I see a couple of issues. Using 4 (or higher) dimensioned arrays probably simplifies your implementation from a mechanics point of view but presents an obstacle when translating to the form PDE Toolbox requires. Second, it is not useful to define the coefficients on a regular grid defined by meshgrid() since PDE Toolbox requires these coefficients at the centroids of each element in the model.
I've included a snippet of MATLAB code below that attempts to address these two issues. For testing, I just used matrices of ones for A1 and A2-- presumably you have the real values for these. You need to take a close look at this code to make sure I haven't made any coding errors and you may need to modify it for your specific case.
c can be defined this way and passed into assempde
A1 = ones(2,2,2,2); A2 = ones(2,2,2,2); % testing only!
c = @(p, t, u, time) calcCmat(p, t, A1, A2, L1);
The function calcCmat is defined by this code:
function c=calcCmat(p, t, A1, A2, L1)
% Triangle point indices
it1=t(1,:);
it2=t(2,:);
it3=t(3,:);
% Find centroids of triangles
X=(p(1,it1)+p(1,it2)+p(1,it3))/3;
eta1 = 0.5*(tanh((X-(L1/3))/0.5)-tanh((X-(2*L1/3))/0.5));
eta2 = 1 - eta1;
% map tensor components to PDE Toolbox c-matrix
ijkl = [1 1 1 1; 1 1 1 2; 1 1 2 2; 1 2 1 1; 1 2 2 1;
2 2 1 1; 1 2 1 2; 1 2 2 2; 2 2 1 2; 2 2 2 2];
c = zeros(10, size(t,2));
for n=1:10
i = ijkl(n,1); j = ijkl(n,2); k = ijkl(n,3); l = ijkl(n,4);
c(n,:) = A1(i,j,k,l)*eta1.^2 + A2(i,j,k,l)*eta2.^2;
end
end
Bill Greene
2013-12-10
>Do you recommend to keep trying to resolve it with this pde tool? Yes, definitely. I'm almost certain that the 10-entry version of the PDE Toolbox c-coefficient will support your case. It is simply a matter of mapping the entries of the 4'th order tensor, C1, to the ten entries in the matrix. I think that the mapping is defined by this vector ijkl = [1 1 1 1; 1 1 1 2; 1 1 2 2; 1 2 1 1; 1 2 2 1; 2 2 1 1; 1 2 1 2; 1 2 2 2; 2 2 1 2; 2 2 2 2]; where each row defines the indices for the 4th order tensor. But I could be wrong since I don't understand the details of your material. If your material is really strange and your c-tensor has no symmetry at all, you can use the 16-row version of the c-matrix which supports a completely general 4th order constitutive tensor in 2D.
Bill
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!