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?

$\endgroup$
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

$\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$ – Professor P. Cosmo Klunk Apr 16 '20 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 Apr 16 '20 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 Apr 17 '20 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 Apr 17 '20 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 Apr 17 '20 at 13:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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