FE Verification using photoelasticity and constrained shape optimization

Prof Veera Sundararaghavan

Dept of Aerospace Engineering

University of Michigan Ann Arbor

Aero 510 Common Project Fiber path optimization

Part 1. Validation of FE model, comparison with photoelastic fringes:

Plexiglass plate Young's modulus E = 3.0 GPa Poisson ratio of 0.3. Plate thickness is 4.98mm. Plane stress y-axis tension of 14.11 MPa

Symmetries (quarter model for part 1)

Ux=0

Uy=0

Ty=14.11 MPa

• I prefer a quad element mesh • Meshes are provided • Use full model for part II and III (fiber path optimization)

Photoelastic experiment When a ray of light passes through a photoelastic material, its electromagnetic wave components are resolved along the two principal stress directions and each component experiences a different refractive index, which is related to the principal stresses.

1 2 where is the phase retardation along ith principal direction 2

i N



   

N can be found from the intensity of analyzer image

Photoelastic fringes

N can also be found from the stress distribution. In Part (i), you will compare the two results.

Experimental image

Finite element mesh provided

• Platehole.plt

• Open in a text editor

• Its in tecplot format

No of nodes

No of elements

coordinates

connectivity

Build upon HW3 code, tecplot is great for plotting

Tecplot format

Add variable names here

Add nodal values here

Rest is same as platehole.plt

Part II

• Fiber path optimization

• Design fiber paths to minimize the peak von Mises stress

• Full model will be simulated with cantilever boundary conditions

Useful book

• Chapter 1 is useful if you have not taken a composites course or Aero 513

• Chapter 1 describes orthotropic stiffness matrix, rotation of stiffness for different fiber angles and plane stress reduction.

• Chapter 6 gives more details on micromechanics if you are interested (not needed for this project, properties will be given)

Fiber properties

• Epoxy+Carbon fiber at 60% fiber volume fraction

• Ef = 241 GPa, vf = 0.2, Em = 3.12 GPa, vm = 0.38

• Micromechanics (chapter 6 of Barbero) will give an orthotropic matrix with equivalent stiffness

• This is for the zero degree fiber

Fiber at an angle (full 3D)

t degrees x1

x2e’2

e’1

l,m,n are cosines of the angles between axes marked

1 ' 'C S

 

cos( ) cos(90 ) cos(90)

cos(90 ) cos( ) cos(90)

cos(90) cos(90) cos(0)

t t

t t





For the above case:

For negative angles (clockwise rotation of fibers), you may use the same formula with t<0. Although overall values of l,m,n will be negative compared to counterclockwise positive convention, T will still be right

Plane stress reduction (full 3D to 2D)

=0 =0 =0

For plane stress the zx,zy and zz stresses are zero

Remove the corresponding rows and columns

Invert to get the D matrix

D matrix

Way to check the composite function (book example)

For a multilayered composite, perform volume averaging: All properties:

Some hints:

NURBs representation

• Non-uniform rational basis spline (NURBS) is a mathematical model commonly used in computer graphics for generating and representing curves and surfaces

• Control points determine the shape of the curve

Design variable: Height of the control points

NURBS surface

Fiber representation on the component

Example These points can move along z axis to generate different fiber topologies

z

Try the function “gnurbs” in the nurbs toolbox on CANVAS site

Function “surfacegen”

• Uses NURBS toolbox, download and add the toolbox to your folder:

https://www.mathworks.com/matlabcentral/file exchange/26390-nurbs-toolbox-by-d-m- spink?focused=5144088&tab=function

Function “surfacegen” provided by instructor

• Input 1: NURBS control point z-locations (5x5 matrix)

• Input 2: Integration points on the specimen (Nintx2 matrix)

• Output: Fiber angles at the integration points (Nintx1 vector of angles in degrees)

Note: Function employs plate with hole geometry of the common project – no need to edit this function.

Example

surfacegen(10*rand(5),[0 20])

An integration point Z coordinate of control points

Output: 68.0946 deg

-10 -5 0 5 10

-30

-20

-10

0

10

20

30

AB 68.0946

Optimization

• Matlab has advanced optimization capabilities

• We will employ the ‘fmincon’ function

• Based on gradient optimization

• MATLAB uses finite differences to find the gradients, the step size can be controlled

Use the following commands to perform optimization:

options = optimoptions('fmincon','FinDiffRelStep',0.1);

X = fmincon(@yourcode,zeros(5,5),[],[],[],[],[],[],[],options)

Your FEM code, should take the 5x5 control point matrix as input and return the max von Mises stress

Example Function maxstress = yourcode(x)

Initial value of the control points. A 5x5 matrix of zeros should give you straight fibers along the y-direction

Step size to change control points during gradient calculation

V. 2015 and above

Fmincon in parallel

options = optimoptions('fmincon','FinDiffRelStep',0.01, 'UseParallel',true);

X = fmincon(@commonproject,zeros(5),[],[],[],[],[],[],[],options)

Solution will be close to optimum with the first two drops - 25+25 runs of the code plus step size control steps

Common project part 2 problem (i)

• Find the fiber topology that minimizes the maximum von Mises stress during y-axis tension (geometry same as part (i) photoelasticity, cantilevered specimen so use full model, fiber properties given in slide 12)

• What is the max allowable load in the optimized specimen if failure stress is 144 MPa, – find % improvement from a straight fiber oriented along the loading direction

-1 0

-5 0

5 1 0

-3 0

-2 0

-1 0 0

1 0

2 0

3 0

AB6 8 .0

9 4 6

Solution from 2017 paper

FEM result

Experimental result

Common project part 2 problem (ii)

• Shear loading

Find the fiber topology that minimizes the maximum von Mises stress during y-axis tension (numbers and geometry same as before, loading is shear instead)

-1 0

-5 0

5 1 0

-3 0

-2 0

-1 0 0

1 0

2 0

3 0

AB6 8 .0

9 4 6

Shear case

• start with zero degree fiber by changing a line in surfacegen

Replace with

Zero degree fiber: stress field in shear

vonmises

387

359

332

304

276

249

221

193

166

138

111

83

55

28

0

Optimum should strengthen this region by changing fiber angles

Whats the stress reduction that you can achieve?

Peak stress region

Reporting requirements

• Final report no more than 6 pages long

• Can be in powerpoint slide format

Future additions

• Multiple layers, vary thicknesses

• thermal loads, buckling instabilities, shell elements

• Multiple NURBS patches for large components, eg. airplane wings

• Manufacturing constraints: minimal interfiber distance, curvature constraints (avoiding sharp angles and voids)