Definition of the c coefficient for assempde

1 次查看（过去 30 天）
laureta
laureta 2013-12-2
回答： Bill Greene 2013-12-10
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

请先登录，再回答此问题。

采纳的回答

Bill Greene
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:
http://www.mathworks.com/help/pde/ug/scalar-functional-form-and-calling-syntax.html
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 个评论
显示 1更早的评论
Bill Greene
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
laureta
laureta 2013-12-5
Hi Bill, My eta vector's components are defined as: eta1(1,1,1,1,1:n2,1:n1) = 0.5*(tanh((X-(L1/3))/0.5)-tanh((X-(2*L1/3))/0.5)); eta2(1,1,1,1,1:n2,1:n1) = 1 - eta1(1,1,1,1,1:n2,1:n1); where n1 and n2 are the xpts and ypts of the triangles of the mesh, L1 is the length of the domain, and [X,Y]=meshgrid(xpts,ypts). The dimensions are choosen to fit with the following loop: C1(i,j,k,l,:,:)=C1(i,j,k,l,:,:)+(eta1(1,1,1,1,:,:).^2).*A1(i,j,k,l) +(eta2(1,1,1,1,:,:).^2).*A2(i,j,k,l) ; A1 and A2 are scalars. As you can imagine, my stiffness tensor is C1. How can I reduce it to a 10 x nr_of_triangles matrix? Should I redefine anything? Thanks in advance!

请先登录，再进行评论。

更多回答（2 个）

Bill Greene
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
  1 个评论
laureta
laureta 2013-12-10
Thanks again for your advices! I have been trying to introduce your changes in my code. Unfortunately, I need to define all the 16 entries in the matrix because the stiffness matrix c is a sum for all entries. The original code was as follows, with Einstein summation:
for i=1:2
for j=1:2
for k=1:2
for l=1:2
for p=1:2
for q=1:2
for r=1:2
for s=1:2
C1(i,j,k,l,:,:)=C1(i,j,k,l,:,:)+(eta1(1,1,1,1,:,:).^2).*(A1(i,p)*A1(j,q)*A1(k,r)*A1(l,s)*C(p,q,r,s,:,:)) ...
+(eta2(1,1,1,1,:,:).^2).*(A2(i,p)*A2(j,q)*A2(k,r)*A2(l,s)*C(p,q,r,s,:,:)) ;
end
end
end
end
end
end
end
end
where A1, A2 are rotation matrices. Sorry for not writing this sum from the beginning! I see it could have been easier... but thanks anyway! Do you recommend to keep trying to resolve it with this pde tool?

请先登录，再进行评论。

Bill Greene
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

类别

在 Help Center 和 File Exchange 中查找有关 General PDEs 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by