Nonlinear finite elements/Homework 4/Solutions

Problem 1: Axially Loaded Bar with Nonlinear Modulus

Given:

Bar of uniform cross-section under axial body load and compressive axial end load ( $R$ ). Nonlinear material with constitutive relation

σ = E_{0} (1 - ε) ε, ε = \frac{d u}{d x} .

Finite element model uses quadratic elements with shape functions

\begin{matrix} N_{1} (x) & = \frac{2 (x_{2} - x) (x_{3} - x)}{h^{2}} \\ N_{2} (x) & = \frac{4 (x - x_{1}) (x_{3} - x)}{h^{2}} \\ N_{3} (x) & = \frac{2 (x - x_{1}) (x - x_{2})}{h^{2}} \end{matrix}

where $h = x_{3} - x_{1}$ .

Solution

Part 1

Template:Font

A plot of the shape functions for an element is shown in Figure 1(a). The shape functions for a mesh are shown in Figure 1(b). If we add the shape functions we get

\begin{matrix} \sum_{i} N_{i} (x) & = \frac{2 (x_{2} - x) (x_{3} - x)}{h^{2}} + \frac{4 (x - x_{1}) (x_{3} - x)}{h^{2}} + \frac{2 (x - x_{1}) (x - x_{2})}{h^{2}} \\ = \frac{2}{h^{2}} [(x_{2} - x) (x_{3} - x) + (x - x_{1}) (x_{3} - x) + (x - x_{1}) (x_{3} - x) + (x - x_{1}) (x - x_{2})] \\ = \frac{2}{h^{2}} [(x_{2} - x_{1}) (x_{3} - x) + (x - x_{1}) (x_{3} - x_{2})] \\ = (\frac{2}{h^{2}}) (\frac{h}{2}) (x_{3} - x_{1}) = 1 ◻ \end{matrix}

File:NFE HW4Prob1 1.png

Figure 1 (a). Shape functions for an element.

File:NFE HW4Prob1 1a.png

Figure 1 (b). Shape functions for a mesh.

Part 2

Template:Font

If we use an isoparametric form of the equations to simplify the integration of the stiffness matrix and load vector terms, the shape functions become

\begin{matrix} N_{1} (ξ) & = - \frac{1}{2} (1 - ξ) ξ \\ N_{2} (ξ) & = (1 + ξ) (1 - ξ) \\ N_{3} (ξ) & = \frac{1}{2} (1 + ξ) ξ . \end{matrix}

The map between $x$ and $ξ$ is

x (ξ) = x_{1} N_{1} (ξ) + x_{2} N_{2} (ξ) + x_{3} N_{3} (ξ); x_{2} = \frac{x_{1} + x_{3}}{2} .

Therefore, the Jacobian of the deformation is

J = \frac{d x}{d ξ} = \frac{h}{2} .

The terms of the stiffness matrix are given by

\begin{matrix} K_{i j} & = \int_{x_{1}}^{x_{3}} A E_{0} (1 - \sum_{k = 1}^{3} u_{k} \frac{d N_{k} (x)}{d x}) \frac{d N_{i} (x)}{d x} \frac{d N_{j} (x)}{d x} d x \\ = \int_{- 1}^{1} A E_{0} [1 - \sum_{k = 1}^{3} u_{k} (J^{- 1} \frac{d N_{k} (ξ)}{d ξ})] (J^{- 1} \frac{d N_{i} (ξ)}{d ξ}) (J^{- 1} \frac{d N_{j} (ξ)}{d ξ}) J d ξ . \end{matrix}

The element stiffness matrix is

𝐊 = \frac{A E_{0}}{3 h^{2}} [\begin{matrix} 7 h + u_{3} + 15 u_{1} - 16 u_{2} & - 8 (h + 2 u_{1} - 2 u_{2}) & h - u_{3} + u_{1} \\ - 8 (h + 2 u_{1} - 2 u_{2}) & 16 (h - u_{3} + u_{1}) & - 8 (h - 2 u_{3} + 2 u_{2}) \\ h - u_{3} + u_{1} & - 8 (h - 2 u_{3} + 2 u_{2}) & 7 h - 15 u_{3} - u_{1} + 16 u_{2} \end{matrix}]

The terms of the tangent stiffness matrix are given by

T_{i j} = K_{i j} + \sum_{m = 1}^{3} \frac{\partial K_{i m}}{\partial u_{j}} u_{m} .

Therefore, the element tangent stiffness matrix is

𝐓 = [\begin{matrix} 7 h + 2 u_{3} + 30 u_{1} - 32 u_{2} & - 8 (h + 4 u_{1} - 4 u_{2}) & h - 2 u_{3} + 2 u_{1} \\ - 8 (h + 4 u_{1} - 4 u_{2}) & 16 (h - 2 u_{3} + 2 u_{1}) & - 8 (h - 4 u_{3} + 4 u_{2}) \\ h - 2 u_{3} + 2 u_{1} & - 8 (h - 4 u_{3} + 4 u_{2}) & 7 h - 30 u_{3} - 2 u_{1} + 32 u_{2} \end{matrix}]

The terms of the element load vector are given by

\begin{matrix} f_{i} & = \int_{x_{1}}^{x_{3}} a x N_{i} (x) d x + {N_{i} (A E \frac{d u}{d x}) |}_{x_{1}}^{x_{3}} \\ = \int_{- 1}^{1} a (\sum_{k = 1}^{3} x_{k} N_{k} (ξ)) N_{i} (ξ) J d ξ + {N_{i} (ξ) R_{i} |}_{- 1}^{1} . \end{matrix}

Therefore, the element load vector is

𝐟 = \frac{a h}{6} [\begin{matrix} x_{1} \\ 2 (x_{1} + x_{3}) \\ x_{3} \end{matrix}] + [\begin{matrix} - R_{1} \\ 0 \\ R_{3} \end{matrix}]

The stress-strain curve in compression is shown in Figure 2.

File:NFE HW4Prob1SigEps.png

Figure 2. Stress-strain curve in compression.

If we use a mesh containing $2$ equally sized elements, we get the solutions shown in Figure 3.

File:NFE HW4Prob1Disp002Elem.png

Figure 3(a). Displacement

File:NFE HW4Prob1Stress002Elem.png

Figure 3(b). Stress

File:NFE HW4Prob1Strain002Elem.png

Figure 3(a). Strain

File:NFE HW4Prob1Residual002Elem.png

Figure 3(b). Error norm

If we use a mesh containing $10$ equally sized elements, we get the solutions shown in Figure 4.

File:NFE HW4Prob1Disp010Elem.png

Figure 4(a). Displacement

File:NFE HW4Prob1Stress010Elem.png

Figure 4(b). Stress

File:NFE HW4Prob1Strain010Elem.png

Figure 4(a). Strain

File:NFE HW4Prob1Residual010Elem.png

Figure 4(b). Error norm

Observations:

The two-node and three-node elements give the same solution.
The solution using two elements shows that the stresses and strains are interpolated linearly within each element and are non longer constant.
There is no jump in the stress/strain at the nodes.
The solution converges quite rapidly - four iterations. This is due to the shape of the stress-strain curve in compression.

Part 2

Template:Font

If we use a mesh containing $100$ equally sized elements, we get the solutions shown in Figure 5.

File:NFE HW4Prob1Disp100Elem.png

Figure 5(a). Displacement

File:NFE HW4Prob1Stress100Elem.png

Figure 5(b). Stress

File:NFE HW4Prob1Strain100Elem.png

Figure 5(a). Strain

File:NFE HW4Prob1Residual100Elem.png

Figure 5(b). Error norm

Observations:

The solution is approximated quite accurately with 10 elements. Using 100 elements does not lead to any significant improvement in the solution.
The solution converges in four iterations unlike the case where a tensile load was applied.

Part 4

Template:Font

The solution converges at all loads unlike in the tension problem. If the load is increased to $R = - 20.1$ , the elements invert and we get a negative Jacobian. Even though we get a solution, that solution is unphysical.

The Matlab code use to compute the results is shown below.


function HW4Prob1_2

  E0 = 10;
  R = -2;
  a = 1;
  W = 1;
  H = 1;
  A = 1;
  L = 1;
  tol = 1.0e-6;

  nElem = 10;
  h = L/nElem;
  nNode = 2*nElem+1;
  for i=1:nNode
    xnode(i) = (i-1)*h/2;
  end
  for i=1:nElem
    node1 = 2*(i-1)+1;
    elem(i,:) = [node1 node1+1 node1+2];
  end

  u = zeros(nNode,1);

  condition = 1;
  condition_res = 1;
  count = 0;
  while (condition > tol)
    count = count + 1

    K = zeros(nNode);
    T = zeros(nNode);
    f = zeros(nNode,1);

    for i=1:nElem
      node1 = elem(i,1);
      node2 = elem(i,2);
      node3 = elem(i,3);
      Ke = elementStiffness(A, E0, xnode, u, node1, node2, node3);
      Te = elementTangentStiffness(A, E0, xnode, u, node1, node2, node3);
      fe = elementLoad(a, xnode, node1, node2, node3);
      K(node1:node3,node1:node3) = K(node1:node3,node1:node3) + Ke;
      T(node1:node3,node1:node3) = T(node1:node3,node1:node3) + Te;
      f(node1:node3) = f(node1:node3) + fe;
    end

    f(nNode) = f(nNode) + R;

    u0 = 0;
    Kred = K(2:nNode,2:nNode);
    Tred = T(2:nNode,2:nNode);
    fred = f(2:nNode);
    ured = u(2:nNode);

    r = Kred*ured - fred;

    Tinv = inv(Tred);

    u_new = ured - Tinv*r;
    condition = norm(u_new - ured)/norm(u_new)
    u = [u0; u_new];
    iter(count) = count;
    normres(count) = norm(r);
    normu(count) = condition;
  end

  figure;
  p01 = semilogy(iter, normres, 'ro-'); hold on;
  p02 = semilogy(iter, normu, 'bo-');
  set([p01 p02], 'LineWidth', 2);
  xlabel('Iterations','FontName','palatino','FontSize',18);
  ylabel('Error Norm','FontName','palatino','FontSize',18);
  set(gca, 'LineWidth', 2, 'FontName','palatino','FontSize',18);
  legend([p01 p02],'|r|','|u_{n+1}-u_n|/|u_{n+1}|');
  axis square;

  figure;
  p0 = plotDisp(E0, A, L, R, a);
  p = plot(xnode, u, 'ro-', 'LineWidth', 2); hold on;
  xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
  ylabel('u(x)', 'FontName', 'palatino', 'FontSize', 18);
  set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);
  legend([p0 p], 'Linear (Exact) ','Nonlinear (FEM)');
  axis square;

  for i=1:nElem
    node1 = elem(i,1);
    node2 = elem(i,2);
    node3 = elem(i,3);
    [xloc(i,:), eps(i,:), sig(i,:)] = ...
    elementStrainStress(E0, xnode, u, node1, node2, node3);
  end

  figure;
  p0 = plotStress(E0, A, L, R, a);
  for i=1:nElem
    p1 = plot(xloc(i,:), sig(i,:), 'r-o','LineWidth',2); hold on;
  end
  xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
  ylabel('\sigma(x)', 'FontName', 'palatino', 'FontSize', 18);
  set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);
  legend([p0 p1], 'Linear (Exact) ','Nonlinear (FEM)');
  axis square;

  figure;
  p0 = plotStrain(E0, A, L, R, a);
  for i=1:nElem
    p1 = plot(xloc(i,:), eps(i,:), 'r-o','LineWidth',2); hold on;
  end
  xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
  ylabel('\epsilon(x)', 'FontName', 'palatino', 'FontSize', 18);
  set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);
  legend([p0 p1], 'Linear (Exact) ','Nonlinear (FEM)');
  axis square;


function [Ke] = elementStiffness(A, E0, xnode, u, node1, node2, node3)

  x1 = xnode(node1);
  x2 = xnode(node2);
  x3 = xnode(node3);
  h = x3 - x1;
  u1 = u(node1);
  u2 = u(node2);
  u3 = u(node3);
  C1 = A*E0/(3*h^2);
  Ke(1,1) = 7*h + u3 + 15*u1 - 16*u2;
  Ke(1,2) = -8*(h + 2*u1 - 2*u2);
  Ke(1,3) = (h - u3 + u1);
  Ke(2,2) = 16*(h - u3 + u1);
  Ke(2,3) = -8*(h - 2*u3 + 2*u2);
  Ke(3,3) = (7*h - 15*u3 - u1 + 16*u2);
  Ke(2,1) = Ke(1,2);
  Ke(3,1) = Ke(1,3);
  Ke(3,2) = Ke(2,3);
  Ke = C1*Ke;

function [Te] = elementTangentStiffness(A, E0, xnode, u, ...
                                        node1, node2, node3)

  x1 = xnode(node1);
  x2 = xnode(node2);
  x3 = xnode(node3);
  h = x3 - x1;
  u1 = u(node1);
  u2 = u(node2);
  u3 = u(node3);
  C1 = A*E0/(3*h^2);
  Te(1,1) = 7*h + 2*u3 + 30*u1 - 32*u2;
  Te(1,2) = -8*(h + 4*u1 - 4*u2);
  Te(1,3) = (h - 2*u3 + 2*u1);
  Te(2,2) = 16*(h - 2*u3 + 2*u1);
  Te(2,3) = -8*(h - 4*u3 + 4*u2);
  Te(3,3) = (7*h - 30*u3 - 2*u1 + 32*u2);
  Te(2,1) = Te(1,2);
  Te(3,1) = Te(1,3);
  Te(3,2) = Te(2,3);
  Te = C1*Te;

function [fe] = elementLoad(a, xnode, node1, node2, node3)

  x1 = xnode(node1);
  x2 = xnode(node2);
  x3 = xnode(node3);
  h = x3 - x1;
  C1 = a*h/6;
  fe(1,1) = x1;
  fe(2,1) = 2*(x1+x3);
  fe(3,1) = x3;
  fe = C1*fe;

function [x, eps, sig] = elementStrainStress(E0, xnode, u, ...
                                             node1, node2, node3)

  x1 = xnode(node1);
  x2 = xnode(node2);
  x3 = xnode(node3);
  h = x3 - x1;
  u1 = u(node1);
  u2 = u(node2);
  u3 = u(node3);

  ximin = -1;
  ximax = 1;
  nxi = 2;
  dxi = (ximax - ximin)/nxi;
  for i=1:nxi+1
    xi = ximin + (i-1)*dxi;
    J = h/2;
    dN1_dxi = xi - 1/2;
    dN2_dxi = -2*xi;
    dN3_dxi = xi + 1/2;
    B = [dN1_dxi/J dN2_dxi/J dN3_dxi/J];
    u = [u1; u2; u3];
    eps(i) = B*u;
    sig(i) = E0*(1-eps(i))*eps(i);
    F = 1 + eps(i);
    if (F < 0)
      display('Material has inverted!');
    end
    N1 = -0.5*(1-xi)*xi;
    N2 = (1+xi)*(1-xi);
    N3 = 0.5*(1+xi)*xi;
    x(i) = x1*N1 + x2*N2 + x3*N3;
  end

function [p] = plotDisp(E, A, L, R, a)

  dx = 0.01;
  nseg = L/dx;
  for i=1:nseg+1
    x(i) = (i-1)*dx;
    u(i) = 1/(6*A*E)*(-a*x(i)^3 + (6*R + 3*a*L^2)*x(i));
  end
  p = plot(x, u, 'LineWidth', 3); hold on;
  xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
  ylabel('u(x)', 'FontName', 'palatino', 'FontSize', 18);
  set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);

function [p] = plotStrain(E, A, L, R, a)

  dx = 0.01;
  nseg = L/dx;
  for i=1:nseg+1
    x(i) = (i-1)*dx;
    eps(i) = 1/(2*A*E)*(-a*x(i)^2 + (2*R + a*L^2));
  end
  p = plot(x, eps, 'LineWidth', 3); hold on;
  xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
  ylabel('\epsilon(x)', 'FontName', 'palatino', 'FontSize', 18);
  set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);

function [p] = plotStress(E, A, L, R, a)

  dx = 0.01;
  nseg = L/dx;
  for i=1:nseg+1
    x(i) = (i-1)*dx;
    sig(i) = 1/(2*A)*(-a*x(i)^2 + (2*R + a*L^2));
  end
  p = plot(x, sig, 'LineWidth', 3); hold on;
  xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
  ylabel('\sigma(x)', 'FontName', 'palatino', 'FontSize', 18);
  set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);

You could alternatively write your own Maple code or use ANSYS for these calculations depending on your interests.

Problem 2: Nonlinear Thermal Conduction

Consider the model chip system shown in Figure 6.

File:ChipHeat.png

Figure 6. Schematic of a model chip.

Assume that heat is being generated in the silicon carbide region of the chip at the rate of 100 W/cm $2$ . The cooling channels contain a fluid at a constant 250 K and (supposedly) can dissipate heat at the rate of 800 W/cm $2$ . The surroundings are at room temperature (300 K).

The thermal conductivity of the silicon carbide and the aluminum nitride are also given in Figure 6. The thermal conductivity of the epoxy resin is 1.5 W/m-K.

File:NFE HW4 p2 mesh.jpg

Figure 7. Meshed model

Part 1

Template:Font

The following ANSYS input is used for this simulation. Due to symmetry, only half of the model will be meshed Figure 7. A constant temperature of 300K is applied around the edge Epoxy resin. Heat is being generated from areas made up by the Silicon Carbide at the rate of 100 W/cm $2$ . Each line forming the cooling channels has constant temperature of 250K.

/prep7

et,1,55
mp,kxx,1,28.634
mp,kxx,2,84.873
mp,kxx,3,1.5

radius = 0.00005
ewidth = 0.0001
eheight = 0.0001
i=0

k,i+1,0,0
k,i+2,radius,0
k,i+3,ewidth,0
k,i+4,ewidth,eheight
k,i+5,0,eheight
k,i+6,0,radius

local,11,1,0,0
k,i+7,radius,45
k,i+8,radius,-45
csys

k,i+9,0,-radius
k,i+10,0,-eheight
k,i+11,ewidth,-eheight

l,i+2,i+3,2
l,i+3,i+4,2
l,i+4,i+5,2
l,i+5,i+6,2

local,11,1,0,0
l,i+6,i+7,2
l,i+7,i+2,2
l,i+2,i+8,2
l,i+8,i+9,2
csys

l,i+9,i+10,2
l,i+10,i+11,2
l,i+11,i+3,2
l,i+7,i+4,2,2
l,i+8,i+11,2

a,i+2,i+3,i+4,i+7
a,i+7,i+4,i+5,i+6
a,i+2,i+3,i+11,i+8
a,i+8,i+11,i+10,i+9

AGEN,7,all,,,7*ewidth,0,0,,1
asel,s,area,,5,28
ARSYM,X,all, , , ,1,0
asel,s,area,,29,52
AGEN,,all,,,7*ewidth*7,,,,,1
asel,all

*do,i,1,37,7
  blc4,ewidth*i,0,ewidth,eheight
  blc4,ewidth*(i+1),0,ewidth,eheight
  blc4,ewidth*(i+2),0,ewidth,eheight
  blc4,ewidth*(i+3),0,ewidth,eheight
  blc4,ewidth*(i+4),0,ewidth,eheight
  blc4,ewidth*i,-eheight,ewidth,eheight
  blc4,ewidth*(i+1),-eheight,ewidth,eheight
  blc4,ewidth*(i+2),-eheight,ewidth,eheight
  blc4,ewidth*(i+3),-eheight,ewidth,eheight
  blc4,ewidth*(i+4),-eheight,ewidth,eheight
*enddo

blc4,ewidth*43,0,ewidth,eheight
blc4,ewidth*43,-eheight,ewidth,eheight
blc4,ewidth*44,0,ewidth,eheight
blc4,ewidth*44,-eheight,ewidth,eheight

*do,i,1,6
  blc4,(i-1)*0.0001,0.0001,0.0001,0.0001
*enddo

*do,i,0,9
  blc4,0.002+(0.0001*i),0.0001,0.0001,0.0001
*enddo

*do,i,0,14
  blc4,0.0005+(i*0.0001),0.0001,0.0001,0.0001
*enddo

*do,i,0,14
  blc4,0.003+(0.0001*i),0.0001,0.0001,0.0001
*enddo

blc4,0.0044,-0.0001,0.0001,0.0001
blc4,0.0044,0,0.0001,0.0001
blc4,0.0044,0.0001,0.0001,0.0001

*do,i,1,51
  blc4,(i-1)*ewidth,2*eheight,ewidth,0.0003
  blc4,(i-1)*ewidth,-eheight,ewidth,-0.0003
*enddo

*do,i,0,5
  blc4,0.0045+(i*0.0001),eheight,ewidth,0.0001
  blc4,0.0045+(i*0.0001),0,ewidth,0.0001
  blc4,0.0045+(i*0.0001),-eheight,ewidth,0.0001
*enddo

aovlap,all

lesize,all,,,2
aatt,1
asel,s,loc,x,0,0.0045
asel,r,loc,y,-0.0001,0.0001
amesh,all
asel,s,loc,x,0,0.0005
asel,a,loc,x,0.002,0.003
asel,r,loc,y,0.0001,0.0002
amesh,all
asel,all
aatt,2
asel,s,loc,x,0.0005,0.002
asel,a,loc,x,0.003,0.0045
asel,r,loc,y,0.0001,0.0002
amesh,all
asel,all
aatt,3
asel,s,loc,x,0,0.005
asel,r,loc,y,-0.0001,-0.0004
amesh,all
asel,s,loc,x,0,0.005
asel,r,loc,y,0.0002,0.0005
amesh,all
asel,s,loc,x,0.0045,0.005
asel,r,loc,y,-0.0001,0.0002
amesh,all
asel,all

ESEL,S,MAT,,1
/color,elem,12,all
esel,all
ESEL,S,MAT,,2
/color,elem,8,all
esel,all
ESEL,S,MAT,,3
/color,elem,4,all
esel,all

nsel,s,loc,x,0,0
dsym,symm,x
nsel,all
ASEL,S,MAT,,2
bfa,all,hgen,1000000
asel,all
lsel,s,loc,x,0,0.005
lsel,r,loc,y,0.0005
DL,all, ,TEMP,300,1
lsel,s,loc,x,0,0.005
lsel,r,loc,y,-0.0004
DL,all, ,TEMP,300,1
lsel,s,loc,y,-0.0004,0.0005
lsel,r,loc,x,0.005
DL,all, ,TEMP,300,1
lsel,all
LSEL,S,LENGTH,,0.3927E-04
DL,all, ,TEMP,250,1
lsel,all

fini

/solu
solve
fini

/post1
plnsol,temp
plnsol,tfsum
fini

Shown in Figure 8 and 9 below are the contour plots of temperature and the vector sum of heat flux. From Figure 9 we see that the cooling channels do not dissipate heat at the rate suggested in the problem statement.

File:NFE HW4 p2 temp.jpg

Figure 8. Temperature of the linear thermal conductivity model (K).

File:NFE HW4 p2 flux.jpg

Figure 9. Flux of the linear thermal conductivity model (W/m-K).

Part 2

Template:Font

Below are the curve-fits for each material.

Silicon Carbide: $k = 2 (1 0^{- 8}) T^{4} - 2 (1 0^{- 5}) T^{3} + 0.0057 T^{2} - 0.7605 T + 87.873$
Aluminum Nitride: $k = - 1 (1 0^{- 9}) T^{4} + 1 (1 0^{- 7}) T^{3} + 0.0003 T^{2} - 0.0941 T + 28.634$

Part 3

Template:Font

The same code that was used in the previous problem is used here. However, the thermal conductivities are tabulated according to their temperature.

/prep7
et,1,55
MPTEMP,,,,,,,,
MPTEMP,1,20
MPTEMP,2,60
MPTEMP,3,100
MPTEMP,4,140
MPTEMP,5,180
MPTEMP,6,220
MPTEMP,7,260
MPTEMP,8,300

MPDATA,KXX,1,,27
MPDATA,KXX,1,,24
MPDATA,KXX,1,,22
MPDATA,KXX,1,,21
MPDATA,KXX,1,,20
MPDATA,KXX,1,,20
MPDATA,KXX,1,,20
MPDATA,KXX,1,,20

MPDATA,KXX,2,,75
MPDATA,KXX,2,,59
MPDATA,KXX,2,,52
MPDATA,KXX,2,,49
MPDATA,KXX,2,,47
MPDATA,KXX,2,,45
MPDATA,KXX,2,,43
MPDATA,KXX,2,,41

mp,kxx,3,1.5

radius = 0.00005
ewidth = 0.0001
eheight = 0.0001
i=0

k,i+1,0,0
k,i+2,radius,0
k,i+3,ewidth,0
k,i+4,ewidth,eheight
k,i+5,0,eheight
k,i+6,0,radius

local,11,1,0,0
k,i+7,radius,45
k,i+8,radius,-45
csys

k,i+9,0,-radius
k,i+10,0,-eheight
k,i+11,ewidth,-eheight

l,i+2,i+3,2
l,i+3,i+4,2
l,i+4,i+5,2
l,i+5,i+6,2

local,11,1,0,0
l,i+6,i+7,2
l,i+7,i+2,2
l,i+2,i+8,2
l,i+8,i+9,2
csys

l,i+9,i+10,2
l,i+10,i+11,2
l,i+11,i+3,2
l,i+7,i+4,2,2
l,i+8,i+11,2

a,i+2,i+3,i+4,i+7
a,i+7,i+4,i+5,i+6
a,i+2,i+3,i+11,i+8
a,i+8,i+11,i+10,i+9

AGEN,7,all,,,7*ewidth,0,0,,1
asel,s,area,,5,28
ARSYM,X,all, , , ,1,0
asel,s,area,,29,52
AGEN,,all,,,7*ewidth*7,,,,,1
asel,all

*do,i,1,37,7
  blc4,ewidth*i,0,ewidth,eheight
  blc4,ewidth*(i+1),0,ewidth,eheight
  blc4,ewidth*(i+2),0,ewidth,eheight
  blc4,ewidth*(i+3),0,ewidth,eheight
  blc4,ewidth*(i+4),0,ewidth,eheight
  blc4,ewidth*i,-eheight,ewidth,eheight
  blc4,ewidth*(i+1),-eheight,ewidth,eheight
  blc4,ewidth*(i+2),-eheight,ewidth,eheight
  blc4,ewidth*(i+3),-eheight,ewidth,eheight
  blc4,ewidth*(i+4),-eheight,ewidth,eheight
*enddo

blc4,ewidth*43,0,ewidth,eheight
blc4,ewidth*43,-eheight,ewidth,eheight
blc4,ewidth*44,0,ewidth,eheight
blc4,ewidth*44,-eheight,ewidth,eheight

*do,i,1,6
  blc4,(i-1)*0.0001,0.0001,0.0001,0.0001
*enddo

*do,i,0,9
  blc4,0.002+(0.0001*i),0.0001,0.0001,0.0001
*enddo

*do,i,0,14
  blc4,0.0005+(i*0.0001),0.0001,0.0001,0.0001
*enddo

*do,i,0,14
  blc4,0.003+(0.0001*i),0.0001,0.0001,0.0001
*enddo

blc4,0.0044,-0.0001,0.0001,0.0001
blc4,0.0044,0,0.0001,0.0001
blc4,0.0044,0.0001,0.0001,0.0001

*do,i,1,51
  blc4,(i-1)*ewidth,2*eheight,ewidth,0.0003
  blc4,(i-1)*ewidth,-eheight,ewidth,-0.0003
*enddo

*do,i,0,5
  blc4,0.0045+(i*0.0001),eheight,ewidth,0.0001
  blc4,0.0045+(i*0.0001),0,ewidth,0.0001
  blc4,0.0045+(i*0.0001),-eheight,ewidth,0.0001
*enddo

aovlap,all

lesize,all,,,2
aatt,1
asel,s,loc,x,0,0.0045
asel,r,loc,y,-0.0001,0.0001
amesh,all
asel,s,loc,x,0,0.0005
asel,a,loc,x,0.002,0.003
asel,r,loc,y,0.0001,0.0002
amesh,all
asel,all
aatt,2
asel,s,loc,x,0.0005,0.002
asel,a,loc,x,0.003,0.0045
asel,r,loc,y,0.0001,0.0002
amesh,all
asel,all
aatt,3
asel,s,loc,x,0,0.005
asel,r,loc,y,-0.0001,-0.0004
amesh,all
asel,s,loc,x,0,0.005
asel,r,loc,y,0.0002,0.0005
amesh,all
asel,s,loc,x,0.0045,0.005
asel,r,loc,y,-0.0001,0.0002
amesh,all
asel,all

ESEL,S,MAT,,1
/color,elem,12,all
esel,all
ESEL,S,MAT,,2
/color,elem,8,all
esel,all
ESEL,S,MAT,,3
/color,elem,4,all
esel,all

KBC,1
nsel,s,loc,x,0,0
dsym,symm,x
nsel,all
ASEL,S,MAT,,2
bfa,all,hgen,1000000
asel,all
lsel,s,loc,x,0,0.005
lsel,r,loc,y,0.0005
DL,all, ,TEMP,300,1
lsel,s,loc,x,0,0.005
lsel,r,loc,y,-0.0004
DL,all, ,TEMP,300,1
lsel,s,loc,y,-0.0004,0.0005
lsel,r,loc,x,0.005
DL,all, ,TEMP,300,1
lsel,all
LSEL,S,LENGTH,,0.3927E-04
DL,all, ,TEMP,250,1
lsel,all

fini

/solu
solve
fini

/post1
plnsol,temp
plnsol,tfsum
 fini

Shown in Figure 10 and 11 below are the contour plots of temperature and the vector sum of heat flux.

File:NFE HW4 p2 nltemp.jpg

Figure 10. Temperature of the nonlinear thermal conductivity model (K).

File:NFE HW4 p2 nlflux.jpg

Figure 11. Flux of the nonlinear thermal conductivity model (W/m-K).

There are many approaches to create the mesh for this problem. You will have to find the way the work best for you. These two codes listed are just two ways that the meshes can be manipulated.

Template:Subpage navbar

Nonlinear finite elements/Homework 4/Solutions

Contents

Problem 1: Axially Loaded Bar with Nonlinear Modulus

Solution

Part 1

Part 2

Part 2

Part 4

Problem 2: Nonlinear Thermal Conduction

Part 1

Part 2

Part 3

Navigation menu

Nonlinear finite elements/Homework 4/Solutions

Problem 1: Axially Loaded Bar with Nonlinear Modulus

Solution

Part 1

Part 2

Part 2

Part 4

Problem 2: Nonlinear Thermal Conduction

Part 1

Part 2

Part 3

Navigation menu

Search