1
$\begingroup$

Following this question and user21 answered, I changed the BCs of the initial problem. The code following bellow.

<< NumericalDifferentialEquationAnalysis`;
Needs["NDSolve`FEM`"];

G = Rationalize[6.894745 10^9]
E1 = Rationalize[26.25 G]; E2 = Rationalize[1.49 G]; G12 = 
 Rationalize[1.04 G]; nu12 = Rationalize[0.28]; nu21 = (E2*nu12)/E1;

t = Rationalize[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 (\[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];
u0 = Rationalize[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/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}, y == 0], 
   DirichletCondition[{v[x, y] == -u0}, y == b]};
omega = Rectangle[{0, 0}, {a, b}];

mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001];

PDEs = FullSimplify[PDEs];

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

Plotting it the solution for U[x,y], it's possible to see two peaks.

U[x,y] plot

Differentiation of x (D[U[x, y], {x, 1}]) potentiates these peaks as can be seen bellow

D[U[x,y],x

These are not expected results. Would anyone know the reason of these peaks? I've been tried to refine the mesh (~ 10000 quad elements) and these peaks keep happening. Am I losing some important thing here?

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

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, 
            Dirichlet

Condition[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]}

Figure 1 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

Figure 2

Improve this answer
$\endgroup$
12
  • $\begingroup$ Alex, you used gammaD = {DirichletCondition[{v[x, y] == u0, u[x, y] == 0}, y == 0], DirichletCondition[{v[x, y] == -u0, u[x, y] == 0}, y == b]};, instead gammaD = {DirichletCondition[{v[x, y] == u0}, y == 0], DirichletCondition[{v[x, y] == -u0, y == b]};. My plots was created using the last one gammaD, as it's in my question. Did you see? $\endgroup$
    – Ponyboy Curtis
    Commented Apr 16, 2020 at 20:56
  • 1
    $\begingroup$ @DiegoMagela Ok, I used a new gammaD and now can see small pick due to boundary condition DirichletCondition[u[x, y] == 0, x == a/2] $\endgroup$
    – Alex Trounev
    Commented Apr 16, 2020 at 22:43
  • 1
    $\begingroup$ @DiegoMagela Yes, it seems like this BC is wrong in your problem. I don't understand how we can physically manage this BC. $\endgroup$
    – Alex Trounev
    Commented Apr 17, 2020 at 11:27
  • 2
    $\begingroup$ @DiegoMagela For alpha=0 there is symmetrical case, so BC is valid, but for $alpha \ne 0$ this BC is not true. $\endgroup$
    – Alex Trounev
    Commented Apr 17, 2020 at 13:27
  • 1
    $\begingroup$ @user21 It is not first time I am getting this kind of error due to mixing of numerical and symbolic calculations. It is better to make all symbolic calculations separately and use final result with 'NDSolve[]'. $\endgroup$
    – Alex Trounev
    Commented Apr 17, 2020 at 13:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.