Revisions to High peaks of FEM solution

added 252 characters in body

Source Link

edited Apr 16, 2020 at 22:38

48.8k
3
51
115

If we change boundary condition as gammaD = {DirichletCondition[{v[x, y] == u0}, y == 0], DirichletCondition[{v[x, y] == -u0}, y == b]};, then the solution looks as follows

added 62 characters in body

Source Link

edited Apr 16, 2020 at 21:10

Alex Trounev

48.8k
3
51
115

Needs["NDSolve`FEM`"];

G = 6.894745 10^9;

E1 = 26.25 G; E2 = 1.49 G; G12 = 
 1.04 G; nu12 = 0.28; nu21 = (E2*nu12)/E1;

t = 0.0050 .0254;
Son = {{1/E1, -nu12/E1, 0}, {-nu21/E2, 1/E2, 0}, {0, 0, 1/G12}}; Qon =
  Inverse[Son];

Q11 = Qon[[1, 1]]; Q12 = Qon[[1, 2]]; Q22 = Qon[[2, 2]]; Q66 = 
 Qon[[3, 3]];

U1 = (3 Q11 + 3 Q22 + 2 Q12 + 4 Q66)/8; U2 = (Q11 - Q22)/
  2; U3 = (Q11 + Q22 - 2 Q12 - 4 Q66)/
  8; U4 = (Q11 + Q22 + 6 Q12 - 4 Q66)/
  8; U5 = (Q11 + Q22 - 2 Q12 + 4 Q66)/8;

alpha = 45.001 (\[Pi]/180);
a = 1; b = 1; d = a Cos[alpha] + b Sin[alpha];
omega = Rectangle[{0, 0}, {a, b}];
mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001, 
   PrecisionGoal -> 5, AccuracyGoal -> 5];
u0 = 0.01;

angle1 = 10; angle0 = 0;
angles = {{angle0, angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {angle0, angle1}};

num = Dimensions[angles][[1]]; h = num*t; pos = Table[0, num + 1];
pos[[1]] = -h/2;
For[i = 2, i <= num + 1, i++, pos[[i]] = pos[[i - 1]] + t];

\[Xi]A = {0, 0, 0, 0, 0};
\[Xi]B = {0, 0, 0, 0, 0};
\[Xi]D = {0, 0, 0, 0, 0};

For[i = 1, i <= num, i++, T0 = angles[[i, 1]];
  T1 = angles[[i, 2]];
  func[s_] := 
   Simplify@((2.0/d) (T1 - T0) Sqrt[(s - d/2)^2] + T0) (\[Pi]/180);
  theta[x_, y_] := alpha + func[x Cos[alpha] + y Sin[alpha]];
  zA = pos[[i + 1]] - pos[[i]]; zB = pos[[i + 1]]^2 - pos[[i]]^2;
  zD = pos[[i + 1]]^3 - pos[[i]]^3;
  V1 = Cos[2 theta[x, y]]; V2 = Sin[2 theta[x, y]];
  V3 = Cos[4 theta[x, y]]; V4 = Sin[4 theta[x, y]];
  \[Xi]a = {1, V1, V2, V3, V4} zA;
  \[Xi]b = {1, V1, V2, V3, V4} zB;
  \[Xi]d = {1, V1, V2, V3, V4} zD;
  \[Xi]A = \[Xi]A + \[Xi]a;
  \[Xi]B = \[Xi]B + \[Xi]b;
  \[Xi]D = \[Xi]D + \[Xi]d;];

mU = {{U1, U2, 0, U3, 0}, {U4, 0, 0, -U3, 0}, {U1, -U2, 0, U3, 0}, {0,
     0, U2/2, 0, U3}, {0, 0, U2/2, 0, -U3}, {U5, 0, 0, -U3, 0}};

mA = mU.\[Xi]A; mB = (mU.\[Xi]B)/2; mD = (mU.\[Xi]D)/3;

A11[x_, y_] = mA[[1]]; A12[x_, y_] = mA[[2]]; A16[x_, y_] = mA[[4]];
A22[x_, y_] = mA[[3]]; A26[x_, y_] = mA[[5]]; A66[x_, y_] = mA[[6]];
D11[x_, y_] = mD[[1]]; D12[x_, y_] = mD[[2]]; D16[x_, y_] = mD[[4]];
D22[x_, y_] = mD[[3]]; D26[x_, y_] = mD[[5]]; D66[x_, y_] = mD[[6]];

Nx[x_, y_] = 
  A11[x, y] D[u[x, y], {x, 1}] + A12[x, y] D[v[x, y], {y, 1}] + 
   A16[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

Ny[x_, y_] = 
  A12[x, y] D[u[x, y], {x, 1}] + A22[x, y] D[v[x, y], {y, 1}] + 
   A26[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

Nxy[x_, y_] = 
  A16[x, y] D[u[x, y], {x, 1}] + A26[x, y] D[v[x, y], {y, 1}] + 
   A66[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

PDEs = {D[Nx[x, y], {x, 1}] + D[Nxy[x, y], {y, 1}], 
   D[Ny[x, y], {y, 1}] + D[Nxy[x, y], {x, 1}]}; 



gammaD = {DirichletCondition[{v[x, y] == u0, u[x, y] == 0}, y == 0], 
   DirichletCondition[{v[x, y] == -u0, u[x, y] == 0}, y == b]};
        
        
        
        {U, V} = NDSolveValue[{PDEs == {0, 0}, gammaD, 
    DirichletCondition[u[x        Dirichlet

Condition[u[x, y] == 0, x == a/2]}, {u, 
        v}, {x, y} \[Element] mesh];

Needs["NDSolve`FEM`"];

G = 6.894745 10^9;

E1 = 26.25 G; E2 = 1.49 G; G12 = 
 1.04 G; nu12 = 0.28; nu21 = (E2*nu12)/E1;

t = 0.0050 .0254;
Son = {{1/E1, -nu12/E1, 0}, {-nu21/E2, 1/E2, 0}, {0, 0, 1/G12}}; Qon =
  Inverse[Son];

Q11 = Qon[[1, 1]]; Q12 = Qon[[1, 2]]; Q22 = Qon[[2, 2]]; Q66 = 
 Qon[[3, 3]];

U1 = (3 Q11 + 3 Q22 + 2 Q12 + 4 Q66)/8; U2 = (Q11 - Q22)/
  2; U3 = (Q11 + Q22 - 2 Q12 - 4 Q66)/
  8; U4 = (Q11 + Q22 + 6 Q12 - 4 Q66)/
  8; U5 = (Q11 + Q22 - 2 Q12 + 4 Q66)/8;

alpha = 45.001 (\[Pi]/180);
a = 1; b = 1; d = a Cos[alpha] + b Sin[alpha];
omega = Rectangle[{0, 0}, {a, b}];
mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001, 
   PrecisionGoal -> 5, AccuracyGoal -> 5];
u0 = 0.01;

angle1 = 10; angle0 = 0;
angles = {{angle0, angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {angle0, angle1}};

num = Dimensions[angles][[1]]; h = num*t; pos = Table[0, num + 1];
pos[[1]] = -h/2;
For[i = 2, i <= num + 1, i++, pos[[i]] = pos[[i - 1]] + t];

\[Xi]A = {0, 0, 0, 0, 0};
\[Xi]B = {0, 0, 0, 0, 0};
\[Xi]D = {0, 0, 0, 0, 0};

For[i = 1, i <= num, i++, T0 = angles[[i, 1]];
  T1 = angles[[i, 2]];
  func[s_] := 
   Simplify@((2.0/d) (T1 - T0) Sqrt[(s - d/2)^2] + T0) (\[Pi]/180);
  theta[x_, y_] := alpha + func[x Cos[alpha] + y Sin[alpha]];
  zA = pos[[i + 1]] - pos[[i]]; zB = pos[[i + 1]]^2 - pos[[i]]^2;
  zD = pos[[i + 1]]^3 - pos[[i]]^3;
  V1 = Cos[2 theta[x, y]]; V2 = Sin[2 theta[x, y]];
  V3 = Cos[4 theta[x, y]]; V4 = Sin[4 theta[x, y]];
  \[Xi]a = {1, V1, V2, V3, V4} zA;
  \[Xi]b = {1, V1, V2, V3, V4} zB;
  \[Xi]d = {1, V1, V2, V3, V4} zD;
  \[Xi]A = \[Xi]A + \[Xi]a;
  \[Xi]B = \[Xi]B + \[Xi]b;
  \[Xi]D = \[Xi]D + \[Xi]d;];

mU = {{U1, U2, 0, U3, 0}, {U4, 0, 0, -U3, 0}, {U1, -U2, 0, U3, 0}, {0,
     0, U2/2, 0, U3}, {0, 0, U2/2, 0, -U3}, {U5, 0, 0, -U3, 0}};

mA = mU.\[Xi]A; mB = (mU.\[Xi]B)/2; mD = (mU.\[Xi]D)/3;

A11[x_, y_] = mA[[1]]; A12[x_, y_] = mA[[2]]; A16[x_, y_] = mA[[4]];
A22[x_, y_] = mA[[3]]; A26[x_, y_] = mA[[5]]; A66[x_, y_] = mA[[6]];
D11[x_, y_] = mD[[1]]; D12[x_, y_] = mD[[2]]; D16[x_, y_] = mD[[4]];
D22[x_, y_] = mD[[3]]; D26[x_, y_] = mD[[5]]; D66[x_, y_] = mD[[6]];

Nx[x_, y_] = 
  A11[x, y] D[u[x, y], {x, 1}] + A12[x, y] D[v[x, y], {y, 1}] + 
   A16[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

Ny[x_, y_] = 
  A12[x, y] D[u[x, y], {x, 1}] + A22[x, y] D[v[x, y], {y, 1}] + 
   A26[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

Nxy[x_, y_] = 
  A16[x, y] D[u[x, y], {x, 1}] + A26[x, y] D[v[x, y], {y, 1}] + 
   A66[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

PDEs = {D[Nx[x, y], {x, 1}] + D[Nxy[x, y], {y, 1}], 
   D[Ny[x, y], {y, 1}] + D[Nxy[x, y], {x, 1}]};

gammaD = {DirichletCondition[{v[x, y] == u0, u[x, y] == 0}, y == 0], 
   DirichletCondition[{v[x, y] == -u0, u[x, y] == 0}, y == b]};



{U, V} = NDSolveValue[{PDEs == {0, 0}, gammaD, 
    DirichletCondition[u[x, y] == 0, x == a/2]}, {u, 
    v}, {x, y} \[Element] mesh];

Needs["NDSolve`FEM`"];

G = 6.894745 10^9;

E1 = 26.25 G; E2 = 1.49 G; G12 = 
 1.04 G; nu12 = 0.28; nu21 = (E2*nu12)/E1;

t = 0.0050 .0254;
Son = {{1/E1, -nu12/E1, 0}, {-nu21/E2, 1/E2, 0}, {0, 0, 1/G12}}; Qon =
  Inverse[Son];

Q11 = Qon[[1, 1]]; Q12 = Qon[[1, 2]]; Q22 = Qon[[2, 2]]; Q66 = 
 Qon[[3, 3]];

U1 = (3 Q11 + 3 Q22 + 2 Q12 + 4 Q66)/8; U2 = (Q11 - Q22)/
  2; U3 = (Q11 + Q22 - 2 Q12 - 4 Q66)/
  8; U4 = (Q11 + Q22 + 6 Q12 - 4 Q66)/
  8; U5 = (Q11 + Q22 - 2 Q12 + 4 Q66)/8;

alpha = 45.001 (\[Pi]/180);
a = 1; b = 1; d = a Cos[alpha] + b Sin[alpha];
omega = Rectangle[{0, 0}, {a, b}];
mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001, 
   PrecisionGoal -> 5, AccuracyGoal -> 5];
u0 = 0.01;

angle1 = 10; angle0 = 0;
angles = {{angle0, angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {angle0, angle1}};

num = Dimensions[angles][[1]]; h = num*t; pos = Table[0, num + 1];
pos[[1]] = -h/2;
For[i = 2, i <= num + 1, i++, pos[[i]] = pos[[i - 1]] + t];

\[Xi]A = {0, 0, 0, 0, 0};
\[Xi]B = {0, 0, 0, 0, 0};
\[Xi]D = {0, 0, 0, 0, 0};

For[i = 1, i <= num, i++, T0 = angles[[i, 1]];
  T1 = angles[[i, 2]];
  func[s_] := 
   Simplify@((2.0/d) (T1 - T0) Sqrt[(s - d/2)^2] + T0) (\[Pi]/180);
  theta[x_, y_] := alpha + func[x Cos[alpha] + y Sin[alpha]];
  zA = pos[[i + 1]] - pos[[i]]; zB = pos[[i + 1]]^2 - pos[[i]]^2;
  zD = pos[[i + 1]]^3 - pos[[i]]^3;
  V1 = Cos[2 theta[x, y]]; V2 = Sin[2 theta[x, y]];
  V3 = Cos[4 theta[x, y]]; V4 = Sin[4 theta[x, y]];
  \[Xi]a = {1, V1, V2, V3, V4} zA;
  \[Xi]b = {1, V1, V2, V3, V4} zB;
  \[Xi]d = {1, V1, V2, V3, V4} zD;
  \[Xi]A = \[Xi]A + \[Xi]a;
  \[Xi]B = \[Xi]B + \[Xi]b;
  \[Xi]D = \[Xi]D + \[Xi]d;];

mU = {{U1, U2, 0, U3, 0}, {U4, 0, 0, -U3, 0}, {U1, -U2, 0, U3, 0}, {0,
     0, U2/2, 0, U3}, {0, 0, U2/2, 0, -U3}, {U5, 0, 0, -U3, 0}};

mA = mU.\[Xi]A; mB = (mU.\[Xi]B)/2; mD = (mU.\[Xi]D)/3;

A11[x_, y_] = mA[[1]]; A12[x_, y_] = mA[[2]]; A16[x_, y_] = mA[[4]];
A22[x_, y_] = mA[[3]]; A26[x_, y_] = mA[[5]]; A66[x_, y_] = mA[[6]];
D11[x_, y_] = mD[[1]]; D12[x_, y_] = mD[[2]]; D16[x_, y_] = mD[[4]];
D22[x_, y_] = mD[[3]]; D26[x_, y_] = mD[[5]]; D66[x_, y_] = mD[[6]];

Nx[x_, y_] = 
  A11[x, y] D[u[x, y], {x, 1}] + A12[x, y] D[v[x, y], {y, 1}] + 
   A16[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

Ny[x_, y_] = 
  A12[x, y] D[u[x, y], {x, 1}] + A22[x, y] D[v[x, y], {y, 1}] + 
   A26[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

Nxy[x_, y_] = 
  A16[x, y] D[u[x, y], {x, 1}] + A26[x, y] D[v[x, y], {y, 1}] + 
   A66[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

PDEs = {D[Nx[x, y], {x, 1}] + D[Nxy[x, y], {y, 1}], 
   D[Ny[x, y], {y, 1}] + D[Nxy[x, y], {x, 1}]}; 



gammaD = {DirichletCondition[{v[x, y] == u0, u[x, y] == 0}, y == 0],DirichletCondition[{v[x, y] == -u0, u[x, y] == 0}, y == b]};
        
        
        
        {U, V} = NDSolveValue[{PDEs == {0, 0}, gammaD, 
            Dirichlet

Condition[u[x, y] == 0, x == a/2]}, {u, 
        v}, {x, y} \[Element] mesh];

Source Link

created Apr 16, 2020 at 20:48

Alex Trounev

48.8k
3
51
115

Don't use exact number since it turns symbolic calculation and probably force LinearSolve to get exact solution. Just put for instance alpha = 45.001 (\[Pi]/180); and then you get desirable result (I deleted line of code you repeat twice with omega = Rectangle[{0, 0}, {a, b}]; mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001];)

Needs["NDSolve`FEM`"];

G = 6.894745 10^9;

E1 = 26.25 G; E2 = 1.49 G; G12 = 
 1.04 G; nu12 = 0.28; nu21 = (E2*nu12)/E1;

t = 0.0050 .0254;
Son = {{1/E1, -nu12/E1, 0}, {-nu21/E2, 1/E2, 0}, {0, 0, 1/G12}}; Qon =
  Inverse[Son];

Q11 = Qon[[1, 1]]; Q12 = Qon[[1, 2]]; Q22 = Qon[[2, 2]]; Q66 = 
 Qon[[3, 3]];

U1 = (3 Q11 + 3 Q22 + 2 Q12 + 4 Q66)/8; U2 = (Q11 - Q22)/
  2; U3 = (Q11 + Q22 - 2 Q12 - 4 Q66)/
  8; U4 = (Q11 + Q22 + 6 Q12 - 4 Q66)/
  8; U5 = (Q11 + Q22 - 2 Q12 + 4 Q66)/8;

alpha = 45.001 (\[Pi]/180);
a = 1; b = 1; d = a Cos[alpha] + b Sin[alpha];
omega = Rectangle[{0, 0}, {a, b}];
mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001, 
   PrecisionGoal -> 5, AccuracyGoal -> 5];
u0 = 0.01;

angle1 = 10; angle0 = 0;
angles = {{angle0, angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {angle0, 
    angle1}, {-angle0, -angle1}, {angle0, angle1}};

num = Dimensions[angles][[1]]; h = num*t; pos = Table[0, num + 1];
pos[[1]] = -h/2;
For[i = 2, i <= num + 1, i++, pos[[i]] = pos[[i - 1]] + t];

\[Xi]A = {0, 0, 0, 0, 0};
\[Xi]B = {0, 0, 0, 0, 0};
\[Xi]D = {0, 0, 0, 0, 0};

For[i = 1, i <= num, i++, T0 = angles[[i, 1]];
  T1 = angles[[i, 2]];
  func[s_] := 
   Simplify@((2.0/d) (T1 - T0) Sqrt[(s - d/2)^2] + T0) (\[Pi]/180);
  theta[x_, y_] := alpha + func[x Cos[alpha] + y Sin[alpha]];
  zA = pos[[i + 1]] - pos[[i]]; zB = pos[[i + 1]]^2 - pos[[i]]^2;
  zD = pos[[i + 1]]^3 - pos[[i]]^3;
  V1 = Cos[2 theta[x, y]]; V2 = Sin[2 theta[x, y]];
  V3 = Cos[4 theta[x, y]]; V4 = Sin[4 theta[x, y]];
  \[Xi]a = {1, V1, V2, V3, V4} zA;
  \[Xi]b = {1, V1, V2, V3, V4} zB;
  \[Xi]d = {1, V1, V2, V3, V4} zD;
  \[Xi]A = \[Xi]A + \[Xi]a;
  \[Xi]B = \[Xi]B + \[Xi]b;
  \[Xi]D = \[Xi]D + \[Xi]d;];

mU = {{U1, U2, 0, U3, 0}, {U4, 0, 0, -U3, 0}, {U1, -U2, 0, U3, 0}, {0,
     0, U2/2, 0, U3}, {0, 0, U2/2, 0, -U3}, {U5, 0, 0, -U3, 0}};

mA = mU.\[Xi]A; mB = (mU.\[Xi]B)/2; mD = (mU.\[Xi]D)/3;

A11[x_, y_] = mA[[1]]; A12[x_, y_] = mA[[2]]; A16[x_, y_] = mA[[4]];
A22[x_, y_] = mA[[3]]; A26[x_, y_] = mA[[5]]; A66[x_, y_] = mA[[6]];
D11[x_, y_] = mD[[1]]; D12[x_, y_] = mD[[2]]; D16[x_, y_] = mD[[4]];
D22[x_, y_] = mD[[3]]; D26[x_, y_] = mD[[5]]; D66[x_, y_] = mD[[6]];

Nx[x_, y_] = 
  A11[x, y] D[u[x, y], {x, 1}] + A12[x, y] D[v[x, y], {y, 1}] + 
   A16[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

Ny[x_, y_] = 
  A12[x, y] D[u[x, y], {x, 1}] + A22[x, y] D[v[x, y], {y, 1}] + 
   A26[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

Nxy[x_, y_] = 
  A16[x, y] D[u[x, y], {x, 1}] + A26[x, y] D[v[x, y], {y, 1}] + 
   A66[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

PDEs = {D[Nx[x, y], {x, 1}] + D[Nxy[x, y], {y, 1}], 
   D[Ny[x, y], {y, 1}] + D[Nxy[x, y], {x, 1}]};

gammaD = {DirichletCondition[{v[x, y] == u0, u[x, y] == 0}, y == 0], 
   DirichletCondition[{v[x, y] == -u0, u[x, y] == 0}, y == b]};



{U, V} = NDSolveValue[{PDEs == {0, 0}, gammaD, 
    DirichletCondition[u[x, y] == 0, x == a/2]}, {u, 
    v}, {x, y} \[Element] mesh];

Visualisation

{DensityPlot[U[x, y], {x, y} \[Element] mesh, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic], 
 DensityPlot[V[x, y], {x, y} \[Element] mesh, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic]}

Stack Exchange Network

Return to Answer