3D elasticity problem with different Young's Modulus in the field

General Discussion
1

Dear FreeFem++ developers

I’m really new to FreeFem++

I want to proceed with simulation in which young’s modulus transforms other objects by dividing them into three parts as shown in the picture below, how should I perform?

I tried to apply modulus differently depending on the label by dividing the label, but I failed.
I hope the block on the left has a high modulus so that it doesn’t deform even if it’s strong and stays in its original shape.
I want that middle block can only deform.

1-1
1-11141×698 16.9 KB

1-2
1-21141×698 18 KB

Below is the code I wrote.

load "msh3"

// Mesh Generation

real S1 = 1;
real H1 = 1;

int move = 1;
int fixed = 10;

border a1(t=0., 1.){x=0; y=(H1-(2*H1)*t)*1.e-3; label = 1;}
border a2(t=0., 1.){x=(S1*t)*1.e-3; y=(-H1)*1.e-3; label = 1;}
border a3(t=0., 1.){x=(S1)*1.e-3; y=(-H1-(-2*H1)*t)*1.e-3; label = 1;}
border a4(t=0., 1.){x=(S1+(-S1)*t)*1.e-3; y=(H1)*1.e-3; label = 1;}

border b1(t=0., 1.){x=S1*1.e-3; y=(H1-(2*H1)*t)*1.e-3; label = 1;}
border b2(t=0., 1.){x=(S1+(S1*t))*1.e-3; y=(-H1)*1.e-3; label = 1;}
border b3(t=0., 1.){x=(2*S1)*1.e-3; y=(-H1-(-2*H1)*t)*1.e-3; label = 1;}
border b4(t=0., 1.){x=(2*S1+(-S1)*t)*1.e-3; y=(H1)*1.e-3; label = 1;}

border c1(t=0., 1.){x=2*S1*1.e-3; y=(H1-(2*H1)*t)*1.e-3; label = 1;}
border c2(t=0., 1.){x=(2*S1+(S1*t))*1.e-3; y=(-H1)*1.e-3; label = 1;}
border c3(t=0., 1.){x=(3*S1)*1.e-3; y=(-H1-(-2*H1)*t)*1.e-3; label = 1;}
border c4(t=0., 1.){x=(3*S1+(-S1)*t)*1.e-3; y=(H1)*1.e-3; label = 1;}

mesh ThS1 = buildmesh (a1(10) + a2(5) + a3(10) + a4(5));
mesh ThS2 = buildmesh (b1(10) + b2(5) + b3(10) + b4(5));
mesh ThS3 = buildmesh (c1(10) + c2(5) + c3(10) + c4(5));

func zmin = -0.1*1.e-3;
func zmax = 0.1*1.e-3;
int MaxLayer = 1;

int[int] rmid = [1, 1];
int[int] rmid2 = [1, 20];
int[int] rmidfix = [1, 10];

int[int] rup = [0, 1];
int[int] rup2 = [0, 20];
int[int] rupfix = [0, 10];

int[int] rdown = [0, 1];
int[int] rdown2 = [0, 20];
int[int] rdownfix = [0, 10];

mesh3 Th1 = buildlayers(ThS1, MaxLayer, zbound=[zmin, zmax], labelmid=rmid, labelup=rupfix, labeldown=rdown);
mesh3 Th2 = buildlayers(ThS2, MaxLayer, zbound=[zmin, zmax], labelmid=rmid, labelup=rup, labeldown=rdown);
mesh3 Th3 = buildlayers(ThS3, MaxLayer, zbound=[zmin, zmax], labelmid=rmid2, labelup=rup2, labeldown=rdown2);

mesh3 Th = Th1 + Th2 + Th3;

plot (Th, wait=1);

// Define Values

real Fx;
real Fy;
real Fz = 200.;

real ym1 = 10.; 		        // Young's Modulus 1
real pr1 = 0.; 				    // Poisson ratio 1
real ym2 = 1.; 		            // Young's Modulus 2
real pr2 = 0.;                  // Poisson ratio 2 

real Mu1 = ym1/2*(1+pr1);
real Lambda1 = ym1*pr1/((1+pr1)*(1-2*pr1));

real Mu2 = ym2/2*(1+pr2);
real Lambda2 = ym2*pr2/((1+pr2)*(1-2*pr2));

// Fespace

fespace Vh(Th,P1);
Vh u1,u2,u3,u4;
Vh v1,v2,v3;
Vh lambda = (label = 1) ? Lambda1 : Lambda2;        // If label = 1 : Low modulus
Vh mu = (label = 1) ? Mu1 : Mu2;                    // If label = 1 : Low modulus

// Macro

real sqrt2=sqrt(2.);
macro epsilon(u1,u2,u3)  [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3))/sqrt2,(dz(u1)+dx(u3))/sqrt2,(dy(u1)+dx(u2))/sqrt2] // EOM
macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM


// Problem

problem Lame([u1,u2,u3],[v1,v2,v3])=
  int3d(Th)(  
	    lambda*div(u1,u2,u3)*div(v1,v2,v3)	
	    +2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) //')
	      )
  - int3d(Th) (
            [Fx, Fy, Fz]' * [v1, v2, v3]
        )
    + on(fixed,u1=0,u2=0,u3=0)
  ;

// Solve

Lame;

u4 = u4 + u3;

Th=movemesh(Th,[x+u1,y+u2,z+u3]);
u1=u1; u2=u2; u3=u3;

plot(Th, prev=1, wait=1 );

Thank you for your cooperation.

Best regards,
Gooyoon

2

Dear Gooyoon,
Maybe you can try to change the logical test in the definition of lambda and mu (the Lamé constants), to use a condition on coordinates instead of label. I am not sure if I understand correctly your code, but S1 and H1 define the three different regions in the mesh? Then one of the logical conditions could be (x<S1).
Furthermore, you may use functions (func) to define the Lamé constants. Now you define them as finite element functions (of type Vh, with P1 FE functions). This choice has implications for the computation of integrals with int3d(Th)(). You can also define a P0 finite element space especially for them, if they are piecewise constants.
Vincent

3

Remark, in FreeFEM the label are only for boundary condition, you can use region name

and use P0 finite element of define lambda and mu

the correction of you script:

Gooyoon.edp (3.2 KB)

4

Dear Vincent

Thank you for your answer.

According to the professor’s answer below, I should define a P0 finite element space especially for the constants (the Lamé constants).
And I think I should use the regions instead of the labels to designate the section.
So I want to designate (region = 1) as low modulus and (region = 10) as high modulus.

Here is the part of the code I defined different constants depend on the regions.

fespace Ph(Th,P0);

Ph lambda = (region = 10) ? Lambda1 : Lambda2;        // If region = 10 : High modulus
Ph mu = (region = 10) ? Mu1 : Mu2;                    // If region = 10 : High modulus

However, the Lamé constants specified in (region = 10) does not apply.
I don’t know why…

Here is the code I ran.

Gooyoon.edp (3.0 KB)

Gooyoon

5

Thank you for your answer, professor!
I’m sorry for the late reply. I’ve tried it in many ways.

I understand I should define different constant values by regions instead of labels.

So I modified the code based on your advice and I changed the situation a little bit.
But modulus setting according to region was not applied.
How do I declare variables (the Lamé constants) according to region?

The result is at below.

1-1
1-1938×559 29.1 KB

1-2
1-2709×476 24.2 KB

This is what I expect.

1-3
1-3849×389 12.2 KB

load "msh3"

// Mesh Generation

real S1 = 1;
real H1 = 1;

int low = 1;
int high = 10;
int fixed = 15;

border a1(t=0., 1.){x=0; y=(H1-(2*H1)*t)*1.e-3; label = fixed;}
border a2(t=0., 1.){x=(S1*t)*1.e-3; y=(-H1)*1.e-3; label = high;}
border a3(t=0., 1.){x=(S1)*1.e-3; y=(-H1-(-2*H1)*t)*1.e-3; label = high;}
border a4(t=0., 1.){x=(S1+(-S1)*t)*1.e-3; y=(H1)*1.e-3; label = high;}

border b1(t=0., 1.){x=S1*1.e-3; y=(H1-(2*H1)*t)*1.e-3; label = high;}
border b2(t=0., 1.){x=(S1+(S1*t))*1.e-3; y=(-H1)*1.e-3; label = low;}
border b3(t=0., 1.){x=(2*S1)*1.e-3; y=(-H1-(-2*H1)*t)*1.e-3; label = high;}
border b4(t=0., 1.){x=(2*S1+(-S1)*t)*1.e-3; y=(H1)*1.e-3; label = low;}

border c1(t=0., 1.){x=2*S1*1.e-3; y=(H1-(2*H1)*t)*1.e-3; label = high;}
border c2(t=0., 1.){x=(2*S1+(S1*t))*1.e-3; y=(-H1)*1.e-3; label = high;}
border c3(t=0., 1.){x=(3*S1)*1.e-3; y=(-H1-(-2*H1)*t)*1.e-3; label = high;}
border c4(t=0., 1.){x=(3*S1+(-S1)*t)*1.e-3; y=(H1)*1.e-3; label = high;}

mesh ThS1 = buildmesh (a1(10) + a2(5) + a3(10) + a4(5),fixedborder=1);    //  to be sure that the point on boundary do not move 
mesh ThS2 = buildmesh (b1(10) + b2(5) + b3(10) + b4(5),fixedborder=1);
mesh ThS3 = buildmesh (c1(10) + c2(5) + c3(10) + c4(5),fixedborder=1);


func zmin = -0.1*1.e-3;
func zmax = 0.1*1.e-3;
int MaxLayer = 1;


int[int] rghigh = [0, 10];
int[int] rglow = [0, 1];

ThS1 = change(ThS1, region=rghigh);
ThS2 = change(ThS2, region=rglow);
ThS3 = change(ThS3, region=rghigh);

mesh3 Th1 = buildlayers(ThS1, MaxLayer, zbound=[zmin, zmax], region=rghigh);
mesh3 Th2 = buildlayers(ThS2, MaxLayer, zbound=[zmin, zmax], region=rglow);
mesh3 Th3 = buildlayers(ThS3, MaxLayer, zbound=[zmin, zmax], region=rghigh);

mesh3 Th = Th1 + Th2 + Th3;

cout << " regions ="<< regions(Th) << endl;
plot (Th, wait=1);

// Define Values

real Fx;
real Fy;
real Fz = 50.;

real ym1 = 10.; 		             // Young's Modulus High
real pr1 = 0.; 				           // Poisson ratio 1
real ym2 = 1.; 		               // Young's Modulus Low
real pr2 = 0.;                   // Poisson ratio 2 

real Mu1 = ym1/2*(1+pr1);
real Lambda1 = ym1*pr1/((1+pr1)*(1-2*pr1));
real Mu2 = ym2/2*(1+pr2);
real Lambda2 = ym2*pr2/((1+pr2)*(1-2*pr2));


// Fespace

fespace Vh(Th,P1);
fespace Ph(Th,P0);
Vh u1,u2,u3,u4;
Vh v1,v2,v3;
Ph lambda = (region = 10) ? Lambda1 : Lambda2;        // If region = 10 : High modulus
Ph mu = (region = 10) ? Mu1 : Mu2;                    // If region = 10 : High modulus

// Macro

real sqrt2=sqrt(2.);
macro epsilon(u1,u2,u3)  [dx(u1),dy(u2),dz(u3),(dz(u2)+dy(u3))/sqrt2,(dz(u1)+dx(u3))/sqrt2,(dy(u1)+dx(u2))/sqrt2] // EOM
macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM

// Problem

problem Lame([u1,u2,u3],[v1,v2,v3])=
  int3d(Th)(  
	    lambda*div(u1,u2,u3)*div(v1,v2,v3)	
	    +2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) ) //')
	      )
  - int3d(Th) (
            [Fx, Fy, Fz]' * [v1, v2, v3]
        )
    + on(fixed,u1=0,u2=0,u3=0)
  ;

// Solve

Lame;

Th=movemesh(Th,[x+u1,y+u2,z+u3]);
u1=u1; u2=u2; u3=u3;

plot(Th, prev=1, wait=1);

Here is the code I modified.

6

Dear Gooyon,
Just a quick answer: maybe you should correct the logical comparison test, the operator is “==” whereas “=” is for assignment.
Please try:

**Ph lambda = (**region == 10) ? Lambda1 : Lambda2; // If region = 10 : High modulus
**Ph mu = (**region == 10) ? Mu1 : Mu2; // If region = 10 : High modulus

Vincent

7

The result with the correct derand a smaller Fz = 50., ici ok

Remember this is linear elasticity so the figure is wrong.

the code :

Gooyoon.edp (3.3 KB)

8

Did you post the right code? The regions don’t appear to reflect the materials
properties until you change as per earlier poster and Fz is still 50 creating
extreme deformation.
It may be easier to see by just dumping to cout,

Ph lambda = (region == 2) ? Lambda1 : Lambda2;        // If label = 1 : Low modulus
Ph mu = (region == 2) ? Mu1 : Mu2;                    // If label = 1 : Low modulus
Vh l1=lambda;
plot(l1,fill=1,wait=1);
// Macro
cout<<" mu "<<mu[]<<endl;
9

I think so , but I can male mistake.

10

All the problems have been resolved! It really worked
Thank you for all your help. I really appreciated.

Best regards,

Gooyoon