Programa Computacional Para Analise e Dimensionamento de Laminados Compositos Em MATLAB

UNIVERSIDADE DE SAO PAULO

ESCOLA DE ENGENHARIA DE SAO CARLOS

DEPARTAMENTO DE ENGENHARIA MECANICA

PROGRAMA COMPUTACIONAL PARA ANALISE E DIMENSIONAMENTO DE

LAMINADOS EM MATERIAIS COMPOSITOS POLIMERICOS REFORCADOS

IMPLEMENTADO EM MATLAB

Gustavo Miranda Guimaraes

Trabalho de Conclusao de Curso apresentado a

Escola de Engenharia de Sao Carlos da

Universidade de Sao Paulo, como requisito para

conclusao do curso de graduacao em Engenharia

Mecanica.

Area de concentracao: Projeto Mecanico

Orientador: Jonas de Carvalho

Sao Carlos, SP

2011

i

Resumo

GUIMARAES, G. M. (2011). Programa computacional para analise de materiais compositos polimericos

reforcados implementado em MATLAB. Trabalho de Conclusao de Curso - Escola de Engenharia de Sao

Carlos, Universidade de Sao Paulo, Sao Carlos, 2011.

Este trabalho tem como objetivo a implementacao de um programa computacional para analise e di-

mensionamento de laminados em materiais compositos polimericos reforcados e sua aplicacao em um

estudo de caso. O programa consiste em uma rotina de calculo baseada na teoria classica dos laminados,

um banco de dados de propriedades de diversos materiais e criterios de falha. Foram implementados

os criterios de falha de maxima tensao, maxima deformacao, Tsai-Hill, Tsai-Wu e Hashin. O programa

foi desenvolvimento em MATLAB. O projeto de uma quilha para um barco a vela e apresentado, para

exemplificar a aplicacao do programa e o uso do acoplamento flexao-torcao.

Palavras chaves: Compositos, Teoria dos Laminados, MATLAB, Quilha.

ii

Abstract

GUIMARAES, G. M. (2011). Programa computacional para analise de materiais compositos polimericos

reforcados implementado em MATLAB. Trabalho de Conclusao de Curso - Escola de Engenharia de Sao

Carlos, Universidade de Sao Paulo, Sao Carlos, 2011.

The objetive of this work is to implement a computer program for design and analysis of fiber rein-

forced composites and use it in a case study. The program consists of a calculation rotine based on

the classical lamination theory, a material property database and failure criteria. The failure criteria

implemented were the maximum stress criterion, maximum strain, Tsai-Hill, Tsai-Wu and Hashin. The

program was developed in MATLAB. A keel blade design is presented, which served as an example for

both the computer program and to show the bending-torsion coupling.

Keywords: Fiber Composites, Laminate Theory, MATLAB, Keel Blade.

iii

Agradecimentos

Em primeiro gostaria de agradecer a famılia por todo o apoio sem o qual nao seria possıvel a conclusao

desta graduacao.

Ao Prof. Dr. Jonas de Carvalho pela orientacao e motivacao em todo o processo de estudo aplicado

ao trabalho.

Ao Prof. Dan Zenkert pelo apoio durante a disciplina e pelos exemplos usados nesse trabalho.

Ao meu amigo Lakshminarayanan Ramamoorthy pela parceria na disciplina.

A Comissao Europeia pela bolsa de estudos concedida para o intercambio pelo programa Erasmus Mundus

External Cooperation Window.

iv

Lista de Figuras

1 Materiais no Boeing 787 Dreamliner [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2 Rotacao de coordenadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

3 Convencao para os carregamentos no laminado. Adaptado de [2] . . . . . . . . . . . . . . 5

4 Deformacoes e tensoes correspondentes com um carregamento Nx. Adaptado de [2] . . . . 5

5 Coordenadas das lamina no laminado. Adaptado de [2] . . . . . . . . . . . . . . . . . . . . 6

6 Deformacoes e tensoes correspondentes com um carregamento Mx. Adaptado de [2] . . . 6

7 Envelope de falha para maxima tensao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

8 Envelope de falha para maxima deformacao . . . . . . . . . . . . . . . . . . . . . . . . . . 10

9 Envelope de falha para Tsai-Hill . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

10 Envelope de falha para Tsai Wu com diferentes valores de F ∗12 . . . . . . . . . . . . . . . . 13

11 Envelope de falha para Hashin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

12 Algoritmo do programa de calculo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

13 Exemplo de um barco a vela com uma quilha . . . . . . . . . . . . . . . . . . . . . . . . . 19

14 Geometria e carregamentos. Adaptado de [2] . . . . . . . . . . . . . . . . . . . . . . . . . 20

15 Momento fletor e torsor em funcao da coordenada x . . . . . . . . . . . . . . . . . . . . . 21

16 Deflexao da quilha para diversos materiais . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

17 Torcao da quilha para diversos materiais . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

18 Termo de acoplamento torcao-flexao para diversos angulos . . . . . . . . . . . . . . . . . . 22

19 Deflexao da quilha em funcao do angulo do laminado . . . . . . . . . . . . . . . . . . . . . 23

20 Torcao da quilha em funcao do angulo do laminado . . . . . . . . . . . . . . . . . . . . . . 23

21 Deflexao e torcao em funcao da variacao da espessura e angulo da lamina . . . . . . . . . 24

22 Hidrofolio e secao vazada equivalente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

23 Deflexao e torcao em funcao da variacaoo da espessura e angulo para a secao vazada . . . 26

24 Espessura otima em funcao da coordenada x . . . . . . . . . . . . . . . . . . . . . . . . . 27

25 Angulo otimo em funcao da coordenada x . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

26 Exemplo do programa de otimizacao de placas tipo sanduıche . . . . . . . . . . . . . . . . 28

v

Lista de Tabelas

1 Possıveis modos de falha para diversos carregamentos. Adaptado de [2] . . . . . . . . . . 8

2 Propriedades mecanicas para o composito carbono (AS4/3501-6) [4] . . . . . . . . . . . . 8

3 Tesoes e deformacoes na juncao quilha-casco para θ = 30 . . . . . . . . . . . . . . . . . . . 24

4 Tesoes e deformacoes para X=1170 e θ = 30 . . . . . . . . . . . . . . . . . . . . . . . . . . 25

vi

Sumario

Resumo i

Abstract ii

Agradecimento iii

1 Introducao e Objetivos 1

1.1 Introducao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2 Revisao Bibliografica 2

2.1 Analise de Tensoes em Laminados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2.2 Criterios de Falha . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3 Metodologia 15

4 Projeto da Quilha 19

4.1 Descricao do problema . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

4.2 Calculo de flexao e torcao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

4.3 Selecao do material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

4.4 Projeto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

4.5 Sugestoes de melhorias no projeto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

5 Conclusoes e Sugestoes 29

A Programa Principal Main.m 31

B Banco de Dados data.m 38

C Funcao para o calculo da matriz de rigidez local e global q.m 39

D Funcao para o calculo da matriz de rigidez local ql.m 40

E Funcao para o calculo da matriz de transformacao T.m 41

1

1 Introducao e Objetivos

1.1 Introducao

Compositos tem sido introduzidos em diversas aplicacoes com inumeras vantagens como reducao de peso,

menor probabilidade por fadiga e corrosao e um potencial maior de integracao estrutural. A industria

aeronautica tem sido pioneira na utilizacao desses materiais. Projetos recentes tem mostrado como a utili-

zacao de materiais compositos. Dois exemplo sao as aeronaves Boeing 787 Dreamliner e Airbus A380. No

Dreamliner, o uso de material composito chega a 50% do peso total, um aumento significativo em relacao

ao Boeing 777 que possui 9% do peso em compositos. A figura 1 mostra a distribuicao de materiais no

Dreamliner.

Figura 1: Materiais no Boeing 787 Dreamliner [1]

Por mais que os compositos apresentam diversas vantagens, seu uso ainda nao e generalizado por di-

versos motivos [2], entre eles a falta de conhecimento de projeto, falta de dados confiaveis, propriedades

altamente dependentes do processo de manufatura, comportamento anisotropico.

1.2 Objetivos

O objetivo desse trabalho e a apresentacao de um programa computacional para analise e dimensiona-

mento de laminados em materiais compositos polimericos reforcados. O calculo de tensoes e baseado na

teoria classica dos laminados de Kirchhoff. Os criterios de falha utilizados foram os de maxima tensao,

maxima deformacao, Tsai-Hill, Tsai-Wu e Hashin. O projeto de uma quilha para um barco a vela e

apresentado, para exemplificar a aplicacao do programa e o uso do acoplamento flexao-torcao.

2

2 Revisao Bibliografica

2.1 Analise de Tensoes em Laminados

2.1.1 Lei de Hooke generalizada

A lei de Hooke, de forma generalizada, pode ser expressa pelo seguinte equacionamento [2]:

σij = Cijklεkl ou ε = Sijklσkl i,j,k,l=1,2,3 (2.1)

Onde C e o tensor de rigidez e S o tensor de flexibilidade. Em notacao vetorial, os tensores de tensao e

deformacao podem ser escritos como:

σ = [σ11, σ22, σ33, σ23, σ31, σ12]t ≡ [σ1, σ2, σ3, σ4, σ5, σ6]t ≡ [σ1, σ2, σ3, τ23, τ31, τ12]t

ε = [ε11, ε22, ε33, 2ε23, 2ε31, 2ε12]t ≡ [ε1, ε2, ε3, ε4, ε5, ε6]t ≡ [ε1, ε2, ε3, γ23, γ31, γ12]t (2.2)

De forma que a lei de Hooke generalizada pode ser expressa por

σ = Cε ou ε = Sσ (2.3)

Onde C e a matriz de rigidez e S a matriz de flexibilidade. A forma generalizada de C e dada por:

C =

C11 C12 C13 C14 C15 C16

C21 C22 C23 C24 C25 C26

C31 C32 C33 C34 C35 C36

C41 C42 C43 C44 C45 C46

C51 C52 C53 C54 C55 C56

C61 C62 C63 C64 C65 C66

(2.4)

As tensoes, deformacoes, matrizes de rigidez e flexibilidade estao em coordenadas locais (1,2,3 ou x′1, x′2, x′3).

Essas coordenadas podem e sao frequentemente mudadas para coordenadas globais (x,y,z ou x1, x2, x3).

No caso de compositos, essa mudanca de coordenada e feita como rotacao no eixo z ou x3, conforme a

figura 2 e a equacao 2.5.

Figura 2: Rotacao de coordenadas

x′1

x′2

x′3

=

cos θ sin θ 0

− sin θ cos θ 0

0 0 1

x1

x2

x3

=

c s 0

-s c 0

0 0 1

x1

x2

x3

(2.5)

3

A mudanca de coordenadas das deformacoes como rotacao no eixo eixo z ou x3 e dada por:

ε11

ε22

ε33

γ23

γ13

γ12

=

c2 s2 0 0 0 sc

s2 c2 0 0 0 −sc

0 0 1 0 0 0

0 0 0 c −s 0

0 0 0 s c 0

−2sc 2sc 0 0 0 c2 − s2

εxx

εyy

εzz

γyz

γxz

γxy

ou εl = T tε (2.6)

No caso da mudanca de coordenadas das tensoes, temos:

σ11

σ22

σ33

τ23

τ13

τ12

=

c2 s2 0 0 0 2sc

s2 c2 0 0 0 −2sc

0 0 1 0 0 0

0 0 0 c −s 0

0 0 0 s c 0

−2sc sc 0 0 0 c2 − s2

σxx

σyy

σzz

τyz

τxz

τxy

ou σl = T−1σ (2.7)

Unindo as equacoes 2.6, 2.7 e 2.3, podemos equacionar a rotacao da matriz de rigidez:

Cl = T−1C(T t)−1

ou C = TClTt

Sl = T tST ou S = (T−1)tSlT−1 (2.8)

Um material anisotropico possui 21 coeficientes independentes do 36 na matriz de rigidez C. Isso

se deve aos dois pontos abaixo:

� O tensores de tensao e deformacao sao simetricos, isto e, σij = σji e εkl = εlk. Isso implica que o

tensor de rigidez e da seguinte forma: Cijkl = Cjikl e Cijkl = Cijlk.

� Se o material for elastico, existe uma funcao de energia de deformacao W tal que

dW = σijεij = Cjiklεkldεij

Que implica qued2W

dεijdεkl= Cijkl e indiferente em relacao a ordem de diferenciacao⇒ Cijkl = Cklij

Caso o material tenha um plano de simetria, o numero de coeficientes independentes na matriz de rigi-

dez C e 13. Esse tipo de material e chamado de material monoclınico. Um material com tres planos

ortogonais de simetria e chamado de material ortotropico. Neste caso, somente 9 coeficientes de C sao

independentes. A matrix de flexibilidade S para um laminado ortotropico e representado por:

S =

1/E1 −ν21/E2 −ν31/E3 0 0 0

−ν12/E1 1/E2 −ν32/E3 0 0 0

−ν13/E1 −ν23/E2 1/E3 0 0 0

0 0 0 1/G23 0 0

0 0 0 0 1/G31 0

0 0 0 0 0 1/G12

(2.9)

4

Se, dada a rotacao por um angulo arbitrario em um dos eixos(ex:eixo 3), as propriedades no plano orto-

gonal (ex:plano 1-2) sao independentes do angulo, esse material e dito transversalmente isotropico, com

5 coeficientes independentes. Por fim, caso o material tenha a mesmas propriedades independentemente

da direcao, esse material e isotropico, caracterizado por dois coeficientes independentes.

2.1.2 Lamina

A lamina e definida como uma fina camada, que no caso de analise de compositos reforcados por fibras,

ortotropica. Dado que a espessura da lamina e bem menor que as outras dimensoes, e possıvel assumir

estado plano de tensoes, isto e, σ3 = τ23 = τ31. Com essas hipoteses, a lei de Hooke se resume a:

ε1

ε2

γ12

=

1/E1 −ν21/E2 0

−ν12/E1 1/E2 0

0 0 1/G12

σ1

σ2

τ12

ou εl = Slσl (2.10)

Que tambem pode ser equacionada por:

σ1

σ2

τ12

=1

1− ν12ν21

E1 ν21E2 0

ν12E1 E2 0

0 0 G12

ε1

ε2

γ12

ou σl = Qlεl (2.11)

A matriz Ql e chamada a matriz de rigidez da lamina nas coordenadas locais. Da mesma forma que as

tensoes, deformacoes e matrizes de rigidez e flexibilidade foram transformadas nas equacoes 2.6, 2.7 e 2.8,

as matrizes Ql, Sl, tensoes e deformacoes tambem podem ser transformadas, por exemplo:

σx

σy

τxy

=

c2 s2 2sc

s2 c2 −2sc

−2sc sc c2 − s2

σ1

σ2

τ12

ou σ = Tσl

2.1.3 Laminado

Um laminado e um empilhamento de laminas formando uma placa. A diferente da lamina, o laminado

possui uma espessura finita de tal forma que tenha rigidez a flexao [2]. O laminado sera analisado

utilizando a Teoria Classica dos Laminados, considerando a seguintes hipoteses [2, 3]:

� Todas as laminas sao macroscopicamente homogeneas e ortotropicas

� Considera-se a adeao perfeita entre as laminas, garantindo assim a continuidade de deslocamentos

� O laminado possui espessura pequena comparada com as dimensoes laterais e essta sob estado plano

de tensoes

� Nao ocorre deslizamento relativo entre as laminas

5

� Hipotese de Kirchhoff: Uma linha originalmente perpendicular ao plano medio continua perpendi-

cular apos a deformacao.

Assim, de forma generalizada:

ε = ε0 + zκ sendo ε = [εx0, εy0, γxy0]t

e κ = [κx, κy, κxy]t

(2.12)

Onde ε e a deformacao no plano medio e κ a curvatura. Dada a convencao para os carregamentos na

figura 3, temos que:

N = [Nx, Ny, Nxy]t

M = [Mx, My, Mxy]t

(2.13)

Figura 3: Convencao para os carregamentos no laminado. Adaptado de [2]

Se uma forca normal Nx for aplicada no laminado, a deformacao sera constante mas ja que a

rigidez na direcao global e diferente para cada lamina, a tensao ira variar de forma incremental, assim

como mostra a figura 4.

Figura 4: Deformacoes e tensoes correspondentes com um carregamento Nx. Adaptado de [2]

O plano de referencia sera transferido para o plano medio do laminado, de acordo com a figura 5.

Cada lamina sera localizada entre dua posicoes zi−1 e zi. A lamina superior possui a primeira coordenada

z = −h/2 e a lamina inferior a ultima coordenada z = h/2.

6

Figura 5: Coordenadas das lamina no laminado. Adaptado de [2]

Usando as coordenadas definidas na figura 5, temos que a tensao media na direcao x no laminado

e:

σx =1

h

h/2∫−h/2

σx dz e a forca na direcao x e Nx =

h/2∫−h/2

σx dz

De forma geral:

N =

h/2∫−h/2

σ dz

Para um laminado carregado com um momento fletor, a situacao e similar, exceto que a deformacao varia

linearmente com a coordenada z, assim como na figura 6.

Figura 6: Deformacoes e tensoes correspondentes com um carregamento Mx. Adaptado de [2]

O momento, de forma geral, e dado por:

M =

h/2∫−h/2

σz dz

O laminado e composto de n laminas com espessuras hi. E possıvel equacionar os carregamentos em

funcao da tensao em cada lamina:

N =

n∑i=1

zi∫zi−1

σi dz e M =

n∑i=1

zi∫zi−1

σiz dz (2.14)

7

Onde σi = [σxi σyi τxyi]. Combinando as equacoes 2.11 e 2.14:

N =

n∑i=1

zi∫zi−1

σi dz =

n∑i=1

Qi

zi∫zi−1

εi dz =

n∑i=1

Qi

zi∫zi−1

(ε0 + zκ) dz

M =

n∑i=1

zi∫zi−1

σiz dz =

n∑i=1

Qi

zi∫zi−1

εiz dz =

n∑i=1

Qi

zi∫zi−1

z (ε0 + zκ) dz

Onde Qi e a matriz de rigidez da lamina i transformada para as coordenadas globais (x,y,z). Em forma

matricial: NM

=

n∑i=1

zi∫zi−1

Qi Qiz

Qiz Qiz2

ε0κ

dz =

n∑i=1

Qi

zi∫zi−1

I Iz

Iz Iz2

dz ·

ε0κ

(2.15)

Onde I e uma matriz identidade de ordem 3. Uma simplificacao e a forma mais usual de se represetar 2.15

e: NM

A B

B D

ε0κ

(2.16)

Tal que

[A, B, D] =

n∑i=1

Qi

[(zi − zi−1),

1

2(z2i − z2i−1),

1

3(z3i − z3i−1)

]A matrix A e a matriz de rigidez extensional, B a matriz de acoplamento e D a matriz de rigidez de

flexao. Reescrevendo de forma expandida a equacao 2.16:

Nx

Ny

Nxy

Mx

My

Mxy

=

A11 A12 A16 B11 B12 B16

A12 A22 A26 B12 B22 B26

A16 A26 A66 B16 B26 B66

B11 B12 B16 D11 D12 D16

B12 B22 B26 D12 D22 D26

B16 B26 B66 D16 D26 D66

εx0

εx0

γxy0

κx

κy

κxy

(2.17)

Na equacao 2.17 e possıvel ver efeitos peculiares nos compositos, em que cargas normais N podem gerar

curvaturas κ devido a matriz B, momentos fletores causarem torcoes (D16 e D26), entre outros. Esses

termos sao chamados termos de acoplamento e serao usados a favor do projeto no estudo de caso da

quilha.

2.2 Criterios de Falha

2.2.1 Introducao

Assim como a rigidez, a resistencia de um laminado e altamente dependente da direcao do carregamento

em relacao as coordenadas do material. A diferenca de resistencia em cada direcao pode ser de algumas

ordens de grandeza. Por exemplo, no material da tabela 2, a relacao entre a resistencia a tracao longi-

tudinal e transversal e de aproximadamente 14 vezes. A dificuldade de e calcular a resistencia de uma

lamina, e consequentemente de uma laminado, se deve ao grande numero de modos de falhas possıveis

para materiais compositos. Alguns exemplo desses modos de falhas sao apresentados na tabela 1. Para

8

prever a resistencia da lamina, alguns criterios de falha sao propostos. Esse criterios podem ser divididos

por criterios interativos e nao interativos, isto e, consideram ou nao a interacao entre os componentes de

tensao ou deformacao [3]. A seguir sao apresentados alguns desses criterios.

Tabela 1: Possıveis modos de falha para diversos carregamentos. Adaptado de [2]

Carregamento Possıvel Modo de Falha

σ1t Tracao Longitudinal Fratura fragil das fibras

Desprendimento de fibras(Pullout) [3]

σ1c Compressao Longitudinal Microflambagem

Falha na interface fibra-matriz

σ2t Tracao Transversal Falha por tracao da matriz

Falha na interface fibra-matriz

σ2c Compressao Transversal Falha por compressao da matriz

Cisalhamento da interface

ˆτ12 Cisalhamento Falha por cisalhamento da matriz

Falha por cisalhamento da interface

Tabela 2: Propriedades mecanicas para o composito carbono (AS4/3501-6) [4]

AS4/3501-6

E1[GPa] 142.04

E2[GPa] 10.34

G12[GPa] 7.17

ν12 0.27

σ1t[MPa] 2275

σ1c[MPa] 1441

σ2t[MPa] 57.23

σ2c[MPa] 227.54

τ12[MPa] 68

ε1t 0,0015

ε1t 0,0083

ε1t 0,005

ε1t 0,021

γ12 0,030

9

2.2.2 Maxima Tensao

O criterio de falha de maxima tensao diz que a tensao em cada direcao principal deve ser menor do que

a resistencia nessa mesma direcao, caso contrario a falha ocorrera. Equacionando, temos:

−σ1c < σ1 < σ1t

−σ2c < σ2 < σ2t

|τ12| < τ12 (2.18)

Cada componente de tensao e independente, sendo entao um criterio nao interativo. A figura 7 mostra

um exemplo deste criterio, no plano σ1 − σ2 para o material da tabela 2. Nota-se que o envelope de

falha definido pelo criterio e um paralelogramo [3] cujo centro geometrico nao coincide com a origem. A

extensao na direcao σ1 se deve a resistencia elevada nesta direcao.

Figura 7: Envelope de falha para maxima tensao

2.2.3 Maxima Deformacao

De forma analoga ao criterio de maxima tensao, o criterio de maxima deformacao diz que deformacao em

cada direcao principal deve ser menor do que a deformacao limite nessa mesma direcao, caso contrario a

falha ocorrera. Equacionando, temos:

−ε1c < ε1 < ε1t

−ε2c < ε2 < ε2t

|γ12| < γ12 (2.19)

Assim como o criterio de maxima tensao, o criterio de maxima deformacao tambem e um criterio nao

interativo. A figura 8 mostra o envelope de falha deste criterio no plano σ1 − σ2. Os criterios de ma-

xima tensao e maxima deformacao serao iguais desde que a relacao tensao-deformacao seja linear ate a

10

falha para todos os modos de falha. A relacao de linearidade usualmente nao ocorre e para a maioria dos

casos o criterio de maxima deformacao apresenta uma melhor correlacao com os dados experimentais [3].

Figura 8: Envelope de falha para maxima deformacao

2.2.4 Tsai-Hill

Os criterios anteriores de maxima tensao e deformacao possuem a desvantagem de nao levar em consi-

deracao a interacao entre carregamentos. Os criterios de Tsai-Hill e Tsai-Wu foram desenvolvidos como

uma tentativa para levar essa interacao em conta.

Existem diversos criterios baseado na energia de distorcao, por exemplo o criterio de von Mises, em que

a deformacao plastica inicia em um material isotropico quando a equacao abaixo for satisfeita.

1

σ2s

{1

2(σx − σy)2 +

1

2(σy − σz)2 +

1

2(σz − σx)2 + 3τ2yz + 3τ2xz + 3τ2xf

}= 1

Hill sugeriu que o criterio de von Mises poderia ser generalizado para materiais anisotropicos acrescen-

tando os parametros F, G, H, L, M e N. Substituindo x, y, z pelos eixos principais 1, 2, 3, o criterio de

Hill e dado por:

H(σ1 − σ2)2 + F (σ2 − σ3)2 +G(σ3 − σ1)2 + Lτ223 +Mτ213 +Nτ212 = 1

Que pode ser reescrito como:

(G+H)σ21 + (F +H)σ2

2 + (F +G)σ23 − 2Hσ1σ2

−2Gσ1σ3 − 2Fσ2σ3 + Lτ223 +Mτ213 +Nτ212 = 1 (2.20)

Os parametros F, G, H, L, M e N podem ser relacionados com os parametros de resistencia convencionais

11

aplicando uma serie de carregamentos unidirecionais. Se somente τ12 atua, entao

N =1

τ212

De forma analoga,

M =1

τ213L =

1

τ223

Se somente σ1 atua:

G+H =1

σ21

De forma analoga,

F +H =1

σ22

F +G =1

σ23

Combinando as tres equacoes anteriores, os parametros F, G e H podem ser expressos pelos parametros

de resistencia do material:

2H =1

σ21

+1

σ22

− 1

σ23

2G =1

σ21

+1

σ23

− 1

σ22

2F =1

σ22

+1

σ23

− 1

σ21

Substituindo os parametros F, G, H, L, M e N na equacao 2.20, temos:

σ1σ21

+σ2σ22

+σ3σ23

−(

1

σ21

+1

σ22

− 1

σ23

)σ1σ2 −

(1

σ21

+1

σ23

− 1

σ22

)σ1σ3

−(

1

σ22

+1

σ23

− 1

σ21

)σ2σ3 +

τ223τ223

+τ213τ213

+τ212τ212

= 1 (2.21)

Assumindo a lamina sendo transversalmente isotropica, isto e, σ2 = σ3, temos:

σ1σ21

+σ2σ22

+σ3σ23

− σ1σ2σ21

− σ1σ3σ21

−(

2

σ22

− 1

σ21

)σ2σ3 +

τ223τ223

+τ213τ213

+τ212τ212

= 1 (2.22)

Assumindo um estado plano de tensao na lamina, isto e,σ3 = τ13 = τ23 = 0, temos:

σ1σ21

− σ1σ2σ21

+σ2σ22

+τ212τ212

= 1 (2.23)

Para aplicar o criterio, os valores corretos de σ1 e σ2 devem ser usados [4]:

σ1 =

σ1t, se σ1 > 0,

σ1c, se σ1 < 0.

σ2 =

σ2t, se σ2 > 0,

σ2c, se σ2 < 0.

A figura 9 mostra o envelope de falha do criterio no plano σ1−σ2. A maior limitacao deste criterio e que

nao ha distincao entre . Alem disso, ele nao e invariante em relacao ao sistemas de coordenadas , sendo

necessario transformar as tensoes atuantes para as coordenadas da lamina [3], [4].

2.2.5 Tsai-Wu

O criterio de Tsai-Wu e um dos criterios mais usados, tambem conhecido como criterio de interacao

quadratica. Como o proprio nome diz, e um metodo interativo. A sua vantagem em relacao ao metodo

12

Figura 9: Envelope de falha para Tsai-Hill

de Tsai-Hill e que e possıvel diferenciar esforcos de tracao e compressao. Alem disso, e invariante em

relacao ao sistema de coordenadas. Sua forma generalizada e [2]:

Fijσiσj + Fiσi onde i, j= 1 a 6 (2.24)

No caso de um estado plano de tensao, i,j = 1,2,6 (σ6 = τ12) e sua forma expandida e:

F11σ21 + F22σ

22 + F66τ

212 + 2F12σ1σ2 + 2F16σ1τ12 + 2F26σ2τ12 + F1σ1 + F2σ2 + F6τ12 = 1 (2.25)

Como ja citado, a resistencia do material deve ser independente do sinal do cisalhamento, que implica

que:

F16 = F26 = F6 = 0 (2.26)

Reescrevendo o criterio temos:

F11σ21 + F22σ

22 + F66τ

212 + 2F12σ1σ2 + F1σ1 + F2σ2 = 1 (2.27)

O criterio e composto por seis parametros, sendo quatro quadraticos e dois lineares. Os parametros

lineares sao responsaveis pela diferenciacao em esforcos de tracao e compressao. Cindo dos seis parametros

sao obtidos de testes mecanicos convencionais. Assim como foi feito para o criterio de Tsai-Hill, aplicando

testes de tracao e compressao na direcao principal temos:

F11σ21t + F1σ1t quando σ1 = σ1t

F11σ21c − F1σ1c quando σ1 = −σ1c (2.28)

Que implica em:

F1 =1

σ1t− 1

σ1cF11

1

σ1tσ1c(2.29)

13

Da mesma forma na direcao transversal temos:

F2 =1

σ2t− 1

σ2cF22

1

σ2tσ2c(2.30)

Aplicando um cisalhamento no plano, temos:

F66 =1

τ212(2.31)

O sexto parametro, F12, requer um carregamento biaxial, que e mais complicado do que um teste uniaxial

convencional. E possıvel determinar limites para este parametro, estabelecendo que o envelope de falha

seja uma curva fechada. Usando o parametro F ∗12, definido na equacao 2.32, seus limites sao −1 < F ∗12 < 1.

F ∗12 =F12√F11F22

(2.32)

Como ja foi citado, obter o valor deste parametro e difıcil, as vezes obtendo um valor fora dos limites

apresentados. Normalmente se assume um valor no intervalo −0.5 < F ∗12 < 0 [2]. Alguns autores [5]

propoem o valor de -0.5, ja outros [6] obtiveram que F ∗12 = 0 melhor representou os ensaios de composito

de fibra de boro/epoxy e fibra de vidro/epoxy. A forma final do criterio e dada por:

σ21

σ1tσ1c+

σ22

σ2tσ2c+τ212ˆτ12

+ 2F12σ1σ2 +σ1σ1t− σ1σ1c

+σ2σ2t− σ2σ2c

= 1 (2.33)

A figura 10 apresenta o envelope de falha para o planoσ1 − σ2 para diferentes valores de F ∗12 utilizando

o material da tabela 2.Diferente dos criterios de maxima tensao e maxima deformacao, o criterio de

Tsai-Wu nao determina o modo de falha diretamente.

Figura 10: Envelope de falha para Tsai Wu com diferentes valores de F ∗12

14

2.2.6 Hashin

O criterio de Hashin, diferente dos criterio de Tsai-Hill e Tsai-Wu, e baseado nos modos de falha. Cada

modo de falha possui um determinado criterio. Como mostrado a seguir, este e um metodo interativo

para determinados modos de falha. Pela propria definicao do metodo, e possıvel diferenciar o modo de

falha. Originalmente, o metodo assumiu que a resistencia na direcao das fibras e somente limitada pela

ruptura das fibras, com valores diferentes em tracao e compressao, sendo formulado como:

σ1 = σ1t quando σ1 > 0

σ1 = σ1c quando σ1 < 0 (2.34)

Falhas devidas a tensao transversal e cisalhamento foram formuladas com o mesmo criterio, ja que sao

modos de falha dominados pela matriz. O criterio e formulado como:(σ2σ2t

)2

+

(τ12ˆτ12

)2

= 1 quando σ2 > 0 (2.35)(σ2σ2t

)2

+

(τ12ˆτ12

)2

= 1 quando σ2 < 0 (2.36)

A figura 11 mostra o envelope de falha no plano σ1 − σ2 para τ12 = 0, o que e equivalente ao criterio de

maxima tensao.

Figura 11: Envelope de falha para Hashin

15

3 Metodologia

O programa foi desenvolvido segundo o algoritmo apresentado na figura 12. O programa e divido em 3

arquivos:

� data.m Banco de dados com as propriedades dos materiais

� Main.m Programa principal com as participais funcoes e criterios de falhas

� q.m Funcao para o calculo das matrizes local e global de rigidez da lamina

� ql.m Funcao para o calculo da matriz de rigidez local

� T.m Funcao para o calculo da matriz de transformacao

Figura 12: Algoritmo do programa de calculo

O banco de dados armazena as propriedades para cada material de acordo com a tabela 2. O usuario pode

inserir um novo material em uma nova coluna na matriz do arquivo data.m. As propriedades necessarias

sao:

16

E1 Modulo de elasticidade na direcao 1 σ2c Tensao maxima de compressao na direcao 2

E2 Modulo de elasticidade na direcao 2 τ12 Tensao maxima de cisalhamento

G12 Modulo de cisalhamento no plano 1-2 ε1t Deformacao maxima de tracao na direcao 1

ν12 Coeficiente de Poisson ε1t Deformacao maxima de compressao na direcao 1

σ1t Tensao maxima de tracao na direcao 1 ε1t Deformacao maxima de tracao na direcao 2

σ1c Tensao maxima de compressao na direcao 1 ε1t Deformacao maxima de compressao na direcao 2

σ2t Tensao maxima de tracao na direcao 2 γ12 Cisalhamento maximo

Alem das propriedades do material, o usuario precisa definir qual a sequencia de laminacao, as propri-

edades de cada lamina, os carregamentos e o tipo de criterio de falha. A sequencia de laminacao e as

propriedades de cada lamina sao inseridas no arquivo Main.m da seguinte forma:

1 %% Laminate Properties

2 % Material Angle Thickness

3 layup = [7 26 3.2

4 7 0 1.8

5 7 0 1.8

6 7 26 3.2];

Cada linha da variavel layup corresponde com uma lamina. A primeira entrada e o tipo do material,

que relaciona com o numero da coluna no banco de dados data.m. A segunda entrada e o angulo de

laminacao e a terceira entrada a espessura da lamina em milımetros. O exemplo do codigo corresponderia

ao laminado [26/0]s.

O carregamento e inserido no arquivo Main.m da seguinte forma:

1 %% Loading

2 load =[24.5 %Nx

3 0 %Ny

4 0 %Nxy

5 −7.3 %Mx

6 0 %My

7 −2829.8]; %Mxy

Cada entrada na variavel load corresponde as seguintes variaveis:

Nx Forca normal na direcao x [N/mm] Mx Momento fletor na direcao x [Nmm/mm]

Ny Forca normal na direcao y [N/mm] My Momento fletor na direcao y [Nmm/mm]

Nxy Forca cisalhante no plano x-y [N/mm] Mxy Momento cisalhante x-y [Nmm/mm]

A escolha do tipo de criterio de falha e feito ajustando a variavel criteria. Cada valor corresponde a

um tipo de criterio. Caso o criterio de Tsai-Wu seja usado, e necessario colocar o valor de F∗12 desejado.

A parte do codigo corresponte e:

17

1 %% Type of Failure Criteria

2 criteria=2; % 1 for Max Stress, 2 Max Strain , 3 Tsai−Hill, 4 Tsai−Wu, 5 Hashin

3 f12s=0; % F*12 coefficient for the Tsai−Wu Failure Criteria Recomended values: −0.5<F*12<0

A primeira parte da parte de calculo no arquivo Main.m corresponde a determinar o numero de laminas,

determinar as coordenas z de inıcio e fim de cada lamina e mudar referencia da coordenada z para o plano

medio do laminado. Esses passos sao feitos da seguinte forma:

1 tamanho=size(layup);

2 n=tamanho(1,1); %Total number of laminae

3 for i=1:n %Stacking lamminae

4 mtrl=layup(i,1);

5 teta=layup(i,2);

6 hint=layup(i,3);

7 z1(i+1)=z1(i)+hint;

8 end

9 h=z1(n+1,1); %Total thickness of the laminate

10 z1=z1−0.5*h*ones(n+1,1);%Corrected Stack Coordinates

Em seguida, a matriz ABD e calculada. A matriz de rigidez nas coordenadas globais e calculada na

funcao q.m. Primeiro e calculada a matriz de rigidez nas coordenadas locais. Em seguida, a matriz

de transformacao e calculada, dado o angulo da lamina. Finalmente, a matriz e transformada para as

coordenadas globais. O codigo fonte da funcao q.m esta disponıvel no apendice C. A matriz ABD e

calculada por:

1 %% ABD

2 for m=1:n %Sum of all Qi in each part A, B and D.

3 Q=q(m,layup);

4 A= A+(z1(m+1)−z1(m))*Q; % Matrix A − Extension Stiffness Matrix

5 B= B+ (1/2)*((z1(m+1))^2 − (z1(m))^2)*Q; % Matrix B − Bending/Extension Coupling Matrix

6 D= D+ (1/3)*((z1(m+1))^3 − (z1(m))^3)*Q; % Matrix D − Bending Stiffness Matrix

7 end

8 ABD=[A B

9 B D];

Com a matriz ABD calculada e os carregamentos disponıveis, e possıvel calcular as deformacoes no plano

medio do laminado e as curvaturas. Isso e feito da seguinte forma no Main.m:

1 %% Global Strains and Curvatures

2 epiglobal=ABD\load;

As tensoes e deformacoes locais sao calculadas em seguida. Esse calculo e feito em dois pontos para cada

lamina, na coordenada inferior e superior. As funcoes ql.m e T.m sao chamadas nesse ponto, disponıveis

nos apendices D e E, respectivamente. A parte do codigo correspondente para o calculo da tensoe e

18

deformacoes locais e:

1 %% Local Strain and Stresses

2 for i=1:n %% epi=epi0 +kappa Z, 0 for "up", 1 for "down"

3 epil0(1:3,i)=T(i,layup)'*(epiglobal(1:3)+epiglobal(4:6)*z1(i));

4 sigmal0(1:3,i)=ql(i,layup)*epil0(1:3,i);

5 epil1(1:3,i)=T(i,layup)'*(epiglobal(1:3)+epiglobal(4:6)*z1(i+1));


7 end

Com as tensoes e deformacoes locais calculadas, e possıvel aplicar o criterio de falha escolhido, que pode

ser:

� Maxima Tensao

� Maxima Deformacao

� Tsai-Hill

� Tsai-Wu

� Hashin

O codigo fonte para cada um dos criterios de falha se encontra dentro programa Main.m no apendice A.

19

4 Projeto da Quilha

4.1 Descricao do problema

Como exemplo de calculo utilizando o programa, foi desenvolvido um projeto de uma quilha de um barco

a vela de corrida, de acordo com exemplo disponıvel em [2]. A quilha consiste em um bulbo de chumbo

preso na ponta de uma lamina. O bulbo serve de contra-peso, estabilizando o barco quando este se inclina

devido as esforcos na vela. A quilha tambem funciona como um hidrofolio, gerando sustentacao via efeito

hidrodinamico. Quanto mais longa a quilha, mais eficiente ela e, tanto para o efeito de contra peso quanto

para a sustentacao.

Figura 13: Exemplo de um barco a vela com uma quilha

A geometria, o carregamento e a requisitos do projeto sao:

L = 4000 mm B = 1300 mm e = 400 mm

q = 0.005 N/mm2 m=3750 kg

Flexao maxima da quilha com angulo de inclinacao de 30°: δ = 150mm

Angulo de torcao da quilha com angulo de inclinacao de 30°: φ = 0.1

Coeficiente de seguranca igual a 3 com angulo de inclinacao de 30°

A figura 14 mostra um modelo simplificado da geometria e carregamentos na quilha. O centro de gravidade

do bulbo esta deslocado uma distancia e da linha de centro da quilha. O peso da lamina da quilha sera

desprezado nos calculos. O esforco hidrodinamico sera aproximado por uma carga uniforme q. Dado

que o CG nao se encontra na linha de centro, a quilha sera submetida tanto a uma flexao quanto uma

torcao. E desejavel que o angulo de torcao da quilha seja minimizado mas ao mesmo tempo utilizando

uma quilha o mais leve quanto possıvel. Uma das vantagens de se utilizar um material composito para

esse tipo de problema e usufruir do acoplamento flexao-torcao que o material pode oferecer. Assume-se

que o caso e muito mais rıgido que a quilha, de forma que a mesma pode ser aproximada por uma viga

20

Figura 14: Geometria e carregamentos. Adaptado de [2]

em balanco. A secao transversal da quilha sera aproximada por uma placa retangular macica, por mais

que a secao correta seria um hidrofolio.

4.2 Calculo de flexao e torcao

Para se calcular a flexao e torcao da quilha, assume-se que o laminado seja simetrico de modo que B=0,

tal que M = Dκ. Apartir da figura 14, e possıvel equacionar os momentos M:

Mx =P (x− L)

b+q(L2 + x2 − 2Lx)

2

My ≈ 0

Mxy =−T2b

(4.1)

Onde P = mg sin θ e T = emg sin θ. As curvaturas κ sao definidas por:

κx = −∂2w

∂x2

κy = −∂2w

∂y2

κxy = −2∂2w

∂x∂y(4.2)

Que podem ser calculadas por κ = D−1M . Integrando as curvaturas, e possıvel obter a deflexao e

torcao da quilha:

δ = w(L, 0) = d11

(−qL4

8+PL3

b3

)+ d16

(TL2

4b

)(4.3)

φponta = d16

(−qL3

6+PL2

2b

)+ d16

(TL

2b

)(4.4)

Os termos d sao os termos da matrix D−1. o termo d16 e responsavel pelo acoplamento flexao-torcao.A

distribuicao do momento fletor e torsor e mostrada na figura 15.

21

Figura 15: Momento fletor e torsor em funcao da coordenada x

4.3 Selecao do material

A selecao de material foi baseada em propriedades mecanicas de sistemas pre-impregnados disponıveis

na literatura [2]. Os materiais foram comparados segundo a deflexao e torcao, para uma lamina com

espessura fixa e angulo variavel, utilizando o mesmo carregamento da quilha de acordo com as figuras 16

e 17. Dada essa comparacao , o material com a melhores propriedades seria o composito com fibras de

boro mas o composito com fibras de carbono HM sera utilizado ja que foi definido no problema que esse

ja era o material do barco. As propriedades desse material sao mostradas na tabela 2.

Figura 16: Deflexao da quilha para diversos materiais

22

Figura 17: Torcao da quilha para diversos materiais

4.4 Projeto

Secao macica

Como ja descrito na secao de calculo de flexao e torcao, o termo d16 e responsavel pelo acoplamento

flexao-torcao. Para o composito com fibras de carbono HM, esse termo e maximo para angulos ±36,

como mostrado na figura 18.

Figura 18: Termo de acoplamento torcao-flexao para diversos angulos

Um laminado inicial foi proposto com somente um angulo, com o carregamento baseado no

angulo de inclinacao da quilha de θ = 30. A variacao da deflexao e da torcao de acordo com o angulo do

laminado sao mostradas nas figuras 19 e 20. A deflexao e mınima para angulos proximos a 0° e maxima

23

para angulos proximos de 90°. Na figura 20 e possıvel ver que para os angulos de 12.2° e 68° a torcao da

quilha e nula. Entre esses dois angulos, a deflexao e menor para 12.2°, assim esse sera o angulo adotado

para o laminado. Para atender o requisito de deflexao maxima de 150mm, e necessario utilizar a espessura

mınima de 43.05mm.

Figura 19: Deflexao da quilha em funcao do angulo do laminado

Figura 20: Torcao da quilha em funcao do angulo do laminado

Dado que as condicoes de carregamento mudam de acordo com o angulo de inclinacao θ e outros

carregamentos nao previstos podem ocorrer, usar um laminado com um so angulo de laminacao nao e

o projeto mais adequado. Um segundo laminado foi proposto, neste caso um laminado simetrico de 4

laminas com o angulo de 12.2° na lamina exterior. O angulo da lamina interna e as espessuras foram

variadas ate atender os requisitos. A figura 21 mostra os requisitos de deflexao e torcao em funcao da

24

variacao da espessura e angulo da lamina. A interseccao entre os requisitos define os valores que os aten-

dem. O laminado utilizado final foi [12.218.45mm/03.05mm]s, com a deflexao de 150mm e torcao de 0.0414°.

Figura 21: Deflexao e torcao em funcao da variacao da espessura e angulo da lamina

Alem dos requisitos de rigidez, a quilha deve resistir aos carregamentos sem falhar. As tensoes e

deformacoes na juncao quilha-casco sao mostradas na tabela 3. O tensoes e deformacoes tambem foram

analisadas para a posicao X=1170, de acordo com a tabela 4.

Tabela 3: Tesoes e deformacoes na juncao quilha-casco para θ = 30

Lamina 1 Lamina 2 Lamina 3 Lamina 4

σ1 55,712 10,19 -8,284 -54,734

σ2 -1,378 -0,258 0,255 1,413

τ12 -2,83 -0,211 0,052 2,644

ε1 0,000310 0,000057 -0,000046 -0,000300

ε1 -0,000220 -0,000041 0,000038 0,000222

γ12 -0,000390 -0,000029 0,000007 0,000369

Utilizando os criterios de falha disponıveis no programa de calculo o laminado nao falha, porem o coefici-

ente de seguranca de acordo com o criterio de maxima deformacao e de somente 1.68, abaixo do requisito

de 3. Um reprojeto do laminado foi feito utilizando [12.225.8mm/03mm]s, que possui um coeficiente de

seguranca de 3. deflexao de 62.52 mm e uma torcao de 0.007°, satisfazendo todos os requisitos.

25

Tabela 4: Tesoes e deformacoes para X=1170 e θ = 30

Lamina 1 Lamina 2 Lamina 3 Lamina 4

σ1 55,712 10,19 -8,284 -54,734

σ2 -1,378 -0,258 0,255 1,413

τ12 -2,83 -0,211 0,052 2,644

ε1 0,000310 0,000057 -0,000046 -0,000300

ε1 -0,000220 -0,000041 0,000038 0,000222

γ12 -0,000390 -0,000029 0,000007 0,000369

Secao Vazada

Uma das simplificacoes apresentadas anteriormente e que a quilha e uma laminado macico. E vantajoso

utilizar um secao vazada, de modo a aumentar a rigidez sem aumentar o peso da estrutura, aproximada-

mente como a uma estrutura do tipo sanduıche. A secao sera aproximada por um laminado de espessura

total de 105 mm, para ter uma secao transversal proxima a do hidrofolio, como apresentado na figura 22.

Figura 22: Hidrofolio e secao vazada equivalente

Assim como no projeto com secao macica, foram utilizadas duas laminas por lado. O lami-

nado inicial [26/0]s utilizado, variando-se a espessura de cada lamina. Pela figura 23, o laminado

[263,2mm/01,8mm]s foi escolhido, com uma secao vazada de 94.6 mm. Esse laminado e 2.9 vezes mais

leve que o laminado macico. a deflexao foi de 50.24mm, torcao de 0.052° e coeficiente de seguranca de

3. Vale comentar que essa e uma aproximacao muito simplificada e que uma estrutura do tipo sanduıche

seria mais elaborada, levando em consideracao a deflexao devido ao cisalhamento e o peso do material

utilizado no nucleo.

26

Figura 23: Deflexao e torcao em funcao da variacaoo da espessura e angulo para a secao vazada

4.5 Sugestoes de melhorias no projeto

Para obter um projeto mais em elaborado, seria possıvel variar tanto a espessura quanto o angulo de

laminacao em funcao da coordenada x. Como mostrado na figura 15, o momento fletor varia em funcao

de x, de forma de para cada posicao existira um laminado otimo. Somente a respeito de Mx, a espessura

do laminado deveria variar de acordo com a figura 24.

27

Figura 24: Espessura otima em funcao da coordenada x

E possıvel tambem variar o angulo de laminacao de forma a obter um angulo de torcao nulo por

toda a quilha. Em cada coordenada, uma combinacao especifica de Mx e Mxy resulta em um angulo de

laminacao para torcao nula, como mostrado na figura 25. Para posicoes acima de X=3550 nao existe um

angulo de laminacao tal que a torcao seja nula. Isso ocorre porque nessa regiao a torcao e totalmente

dominante e a flexao quase nula, tornando impossıvel zerar a torcao pelo acoplamento flexao-torcao.

Dado que a flexao e dominante proxima a juncao quilha-casco e a torcao proxima a ponta, seria adequado

angulos de laminacao perto de 0° na juncao e 45° na ponta.

Figura 25: Angulo otimo em funcao da coordenada x

Em paralelo com o programa para calculo de compositos, foi desenvolvido um programa para cal-

28

culo de estruturas do tipo sanduıche e um programa para otimizacao com base nas rotinas de otimizacao

fmincon do MATLAB. Um exemplo de otimizacao de placas tipo sanduıche e mostrado na figura 26. E

possıvel adaptar o programa para projetar a quilha, levando em consideracao as deformacoes cisalhantes

no nucleo, restricoes de deflexao e torcao e criterios de resistencia. E possıvel parametrizar a espessura

e angulo de laminacao em funcao da coordenada x. Vale lembrar que a utilizacao dessas rotinas sao de-

pendentes do ponto inicial, entao seria interessante avaliar outras rotinas de otimizacao, como algoritmos

geneticos e recozimento simulado.

Figura 26: Exemplo do programa de otimizacao de placas tipo sanduıche

29

5 Conclusoes e Sugestoes

Foi possıvel equacionar e criar codigos gerais para aplicacao em projeto e analise em compositos. Diversos

criterios de falha foram aplicados, sendo possıvel comparar cada criterio. O projeto da quilha do barco

a vela utilizou o programa e exemplificou o acoplamento flexao-torcao, neste caso usado a favor do projeto.

O codigo poderia ser mais completo se os seguintes topicos estudados fossem adicionados:

� Analise de falha progressiva

� Concentracao de tensao em placas anisotropicas

� Tensoes interlaminares

� Efeitos higrotermicos

A analise progressiva, principalmente associada a um software de elementos finitos, e de importancia para

calcular a progressao da falha no laminado e determinar a sua resistencia final.

Descontinuidades geometricas geram concentracoes de tensoes. Dependendo da sequencia de lamina-

cao e do material, o fator de concentracao de tensao pode ser muito elevado, superior ao caso equivalente

com materiais isotropicos. Esse efeito pode ser calculado analiticamente utilizando as equacoes de con-

centracao de tensao em placas anisotropicas de Lekhnitskii [7].

Tensoes interlaminares sao responsaveis por delaminacoes. Essas tensoes possuem origem diversas, como

carregamentos transversais, diferenca do coeficiente de Poisson entre as laminas, efeitos de borda, furos,

etc. Em alguns casos, e possıvel calcular analiticamente as tensoes interlaminares, levando-as em conta

no projeto do laminado. Para casos mais complexos, uma simulacao numerica e mais adequada.

Os efeitos higrotermicos sao necessarios para calcular tanto tensoes residuais devido ao processo de fa-

bricacao quando efeitos com empenamento devido a mudanca de temperatura e umidade. As tensoes

residuais devem ser levadas em conta para uma melhor estimativa da resistencia do laminado [4].

30

Referencias

[1] GRIFFITHS, B. Boeing sets pace for composite usage in large civil aircraft. High-Performance Com-

posites, May 2005.

[2] ZENKERT, D.; BATTLEY, M. Foundations of Fibre Composites: Royal Institute of Technology,

2003.

[3] SOUZA, G. P. Avaliacao de criterios de falhas de compositos polimericos reforcados aplicados a vigas

sob carregamento de flexao. Dissertacao (Mestrado) — Escola de Engenharia de Sao Carlos, Universi-

dade de Sao Paulo, 2003.

[4] DANIEL, I.; ISHAI, O. Engineering Mechanics of Composite Materials. 2. ed.: Oxford University

Press, USA, 1994.

[5] HYER, M.; WHITE, S. Stress analysis of fiber-reinforced composite materials: WCB McGraw-Hill,

1998.

[6] NARAYANASWAMI, R.; ADELMAN, H. Evaluation of the tensor polynomial and hoffman strength

theories for composite materials. Journal of Composite Materials, SAGE Publications, v. 11, n. 4,

p. 366, 1977.

[7] LEKHNITSKII, S. G. Anisotropic Plates: Science Publishers Inc, 1968.

[8] HASHIN, Z. Failure criteria for unidirectional fiber composites. Journal of Applied Mechanics, v. 47,

p. 329, 1980.

[9] KIM, Y.; DAVALOS, J.; BARBERO, E. Progressive failure analysis of laminated composite beams.

Journal of Composite Materials, SAGE Publications, v. 30, n. 5, p. 536–560, 1996.

[10] SUN, C. et al. Comparative evaluation of failure analysis methods for composite laminates. Washing-

ton, DC: FAA, Office of Aviation Research, 1996.

[11] TSAI, S.; WU, E. A general theory of strength for anisotropic materials. Journal of composite

materials, SAGE Publications, v. 5, n. 1, p. 58, 1971.

31

A Programa Principal Main.m

1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2 % Trabalho de Conclusão de Curso 2011

3 % Engenharia Mcânica− EESC USP

4 % Gustavo Miranda Guimaraes

5 % [email protected]

6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

7

8 clc

9 clear all


11 % Material Angle Thickness

12 layup = [7 26 3.2

13 7 0 1.8

14 7 0 1.8

15 7 26 3.2];

16 %% Loading

17 load =[24.5 %Nx

18 0 %Ny

19 0 %Nxy

20 −7.3 %Mx

21 0 %My

22 −2829.8]; %Mxy

23 %% Type of Failure Criteria

24 criteria=2; % 1 for Max Stress, 2 Max Strain , 3 Tsai−Hill, 4 Tsai−Wu, 5 Hashin

25 f12s=0; % F*12 coefficient for the Tsai−Wu Failure Criteria Recomended values: −0.5<F*12<0

26


28 tamanho=size(layup);

29 n=tamanho(1,1); %Total number of laminae

30 z1=zeros(n+1,1);

31 sf=zeros(n,3); % saftey factor 1st column for 1st direction ,2nd column for 2nd ...

direction of fibre and 3rd column for shear

32 for i=1:n %Stacking lamminae

33 mtrl=layup(i,1);

34 teta=layup(i,2);

35 hint=layup(i,3);

36 z1(i+1)=z1(i)+hint;

37 end

38 h=z1(n+1,1); %Total thickness of the laminate

39 z1=z1−0.5*h*ones(n+1,1); %Corrected Stack Coordinates

40

41 %% ABD

42 A=zeros(3,3);

43 B=zeros(3,3);

44 D=zeros(3,3);

45

32

46 for m=1:n %Sum of all Qi in each part A, B and D.

47 Q=q(m,layup);

48 A= A+(z1(m+1)−z1(m))*Q; % Matrix A − Extension Stiffness Matrix

49 B= B+ (1/2)*((z1(m+1))^2 − (z1(m))^2)*Q; % Matrix B − Bending/Extension Coupling Matrix

50 D= D+ (1/3)*((z1(m+1))^3 − (z1(m))^3)*Q; % Matrix D − Bending Stiffness Matrix

51 end

52 ABD=[A B

53 B D];

54

55 %% Homogenized Elastic Properties

56

57 Qt=A/h;

58 S=inv(Qt);

59 Ex=1/(S(1,1));

60 Ey=1/(S(2,2));

61 Gxy=1/(S(3,3));

62 nuxy=−Ex*S(2,1);

63

64

65 %% Global Strains and Curvatures

66 epiglobal=ABD\load;

67

68 %% Local Strain and Stresses

69 epil0 = zeros(3,n);

70 sigmal0= zeros(3,n);

71 epil1 = zeros(3,n);

72 sigmal1= zeros(3,n);

73

74 for i=1:n %% epi=epi0 +kappa Z, 0 for "up", 1 for "down"

75 epil0(1:3,i)=T(i,layup)'*(epiglobal(1:3)+epiglobal(4:6)*z1(i));


77 epil1(1:3,i)=T(i,layup)'*(epiglobal(1:3)+epiglobal(4:6)*z1(i+1));


79 end

80

81

82 %% Failure Criteria

83

84 if criteria==1 % Max Stress

85 for i=1:n % counts for all plies

86 mtrl=layup(i,1);

87

88 if (sigmal0(1,i)) > 0

89

90 if ( (abs(sigmal0(1,i)) > data(6,mtrl)) ) %sigma in 1st ...

direction check

91 disp(['failed in 1st direction in tension layer ' , num2str(i) ])

92 else

93 sf(i,1)= data(6,mtrl) / abs(sigmal0(1,i));

33

94 end

95 % else 'OK'

96 else

97 if (abs(sigmal0(1,i)) > data(7,mtrl))

98 disp(['failed in 1st direction in compression layer ' , num2str(i) ])

99 else

100 sf(i,1)=data(7,mtrl)/abs(sigmal0(1,i));

101 end

102

103 end

104

105 if sigmal0(2,i) > 0

106 if ( (abs(sigmal0(2,i)) > data(8,mtrl)) ) % sigma in 2nd direction ...

check

107 disp(['failed in 2nd direction in tension layer ',num2str(i)])

108 else


110 end

111

112 else

113

114 if (abs(sigmal0(2,i)) > data(9,mtrl))

115 disp(['failed in 2nd direction in compression layer ',num2str(i)])

116 else


118 end

119

120 end

121

122 % else 'ok'

123

124

125 if (abs(sigmal0(3,i)) > data(10,mtrl)) % sigma in shear

126 disp(['failed in shear layer ', num2str(i)])

127 else


129 %else 'ok'

130 end

131

132 end

133 criteria

134 end

135

136

137 if criteria == 2 % Max Strain

138 for i=1:n % counts for all plies

139 mtrl=layup(i,1);

140 if epil0(1,i) > 0

141

34

142 if ( (abs(epil0(1,i)) > data(11,mtrl)) ) % epi in 1st ...

direction check

143 disp(['failed in 1st direction in tension layer ' , num2str(i)])

144 else

145 sf(i,1)= (data(11,mtrl)/abs(epil0(1,i)));

146 end

147

148 else

149

150 if (abs(epil0(1,i)) > data(12,mtrl))

151 disp(['failed in 1st direction in compression layer ' , num2str(i)] )

152 else

153 sf(i,1)= (data(12,mtrl)/abs(epil0(1,i)));

154 end

155

156 % else 'OK'

157 end

158 if epil0(2,i)> 0

159 if ( (abs(epil0(2,i)) > data(13,mtrl)) ) % epi in 2nd direction check

160 disp(['failed in 2nd direction in tension layer ',num2str(i)])

161 else

162 sf(i,2)=( data(13,mtrl))/(abs(epil0(2,i)));

163

164 end

165 else

166 if (abs(epil0(2,i)) > data(14,mtrl))

167 disp(['failed in 2nd direction in compression layer ',num2str(i)])

168 else

169 sf(i,2)= data(14,mtrl)/abs(epil0(2,i));

170 end

171

172

173 % else 'ok'

174 end

175

176 if (abs(epil0(3,i)) > data(15,mtrl)) % epi in shear

177 disp(['failed in shear layer ' ,num2str(i)])

178 else

179 sf(i,3)=data(15,mtrl)/abs(epil0(3,i));

180 %else 'ok'

181 end

182

183

184 end

185 criteria;

186 sf

187 end

188

189

35

190 if criteria == 3 % Tsai Hill

191

192 test=zeros(1,n);

193 tsai=zeros(1,n);

194 for i=1:n


196

197 if sigmal0(1,i)>0 %Sigma 1 Sign Check

198 s1=data(6,mtrl);

199 else


201 end

202

203 if sigmal0(2,i)>0 %Sigma 2 Sign Check


205 else


207 end

208

209 t12=data(10,mtrl);

210

211

212 f1= sigmal0(1,i)^2 / (s1^2);

213

214 f2= ( (sigmal0(1,i) *sigmal0(2,i) ) /(s1^2));

215 f3=( ((sigmal0(2,i))^2 ) / (s2^2) );

216 f4=(((sigmal0(3,i))^2)/(t12^2));

217

218

219 TH= (f1−f2+f3+f4); %Tsai−Hill Function

220

221 tsai(i)=TH;

222 if TH>1

223 test(i)=1;

224 else

225 test(i)=0;

226 end

227

228 end

229 T2=sum(test);

230 tsai;

231 if T2>0

232 res=('FAIL')

233 else

234 res=('OK')

235 end

236 sfth=1./tsai;% safety factor tsai hill

237 criteria ;

238 end

36

239

240 if criteria ==4 % Tsai WU

241


243 tsaiw=zeros(1,n);

244 for i=1:n


246 s1t=data(6,mtrl);

247 s1c=data(7,mtrl);

248 s2t=data(8,mtrl);

249 s2c=data(9,mtrl);

250 t12h=data(10,mtrl);

251 f12=f12s*sqrt((1)/(s1t*s1c)*(1)/(s2t*s2c));

252 s1=sigmal0(1,i);

253 s2=sigmal0(2,i);

254 t12=sigmal0(3,i);

255 TW=(s1^2)/(s1t*s1c)+(s2^2)/(s2t*s2c)+(t12/t12h)^2+(s1/sit)−(si/sic)+(s2/s2t)−(s2/s2c)+(2*f12*s1*s2); ...

%Tsai−Wu Function

256 tsaiw(i)=TH;

257 if TH>1

258 test(i)=1;

259 else

260 test(i)=0;

261 end

262

263 end

264 T2=sum(test);

265 if T2>0

266 res=('FAIL')

267 else

268 res=('OK')

269 end

270 sfth=1./tsaiw;% safety factor Tsai Wu

271 end

272

273 if criteria == 5 % Hashin

274


276 htest=zeros(1,n);

277 for i=1:n


279

280 if sigmal0(1,i)>0 %Sigma 1 Signal Check


282 else


284 end

285

286 if sigmal0(2,i)>0 %Sigma 2 Signal Check

37


288 else


290 end

291

292 t12=data(10,mtrl);

293

294

295 hashintest(1)=(abs(sigmal0(1,i))/s1);

296 hashintest(2)=((sigmal0(2,i)/s2)^2+(sigmal0(3,i)/t12)^2);

297

298

299 Hashin= max(hashintest);

300

301 htest(i)=TH;

302 if TH>1

303 test(i)=1;

304 else

305 test(i)=0;

306 end

307

308 end

309 T2=sum(test);

310 if T2>0

311 res=('FAIL')

312 else

313 res=('OK')

314 end

315 sfth=1./htest;% safety factor for the Hashin criteria

316 end

38

B Banco de Dados data.m

1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





6 %

7 % Banco de dados de

8 % ZENKERT, D.; BATTLEY, M. Foundations of Fibre Composites.

9 % Royal Institute of Technology, 2003.

10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

11 function [Q]=data(var,mat);

12 %% Material Properties

13

14 % Lamina 1 2 3 4 5 6 7

15 % Fibre E−glass Boron Carbon (HT) Carbon (IM) Kevlar 49 Carbon Carbon (HM)

16

17 data = [0.127 0.220 0.129 0.129 0.132 0.127 13.05 52.3% Thickness [mm]1

18 40000 210000 136000 151000 75000 147000 181000 0.001% E1 [MPa]2

19 9800 20000 10000 9400 6000 9000 10300 0.001% E2 [MPa]3

20 2800 6000 5200 4800 2000 3300 7170 0.001% G12 [MPa]4

21 0.3 0.3 0.3 0.31 0.34 0.31 0.28 0.0001% nu12 5

22 1100 1400 1800 2260 1400 2260 1500 9000% Sigma1t [MPa]6

23 600 2800 1200 1200 280 1200 1500 9000% Sigma1c [MPa]7

24 20 80 40 50 30 50 40 9000% Sigma2t [MPa]8

25 140 280 220 190 140 190 246 9000% Sigma2c [MPa]9

26 70 120 80 100 60 100 68 9000% Tau12c [MPa]10

27 0.028 0.007 0.013 0.015 0.019 0.015 0.0083 1% Epsilon1t 11

28 0.015 0.013 0.009 0.008 0.004 0.008 0.0083 1% Epsilon1c 12

29 0.002 0.004 0.004 0.005 0.005 0.005 0.0039 1% Epsilon2t 13

30 0.014 0.014 0.021 0.020 0.023 0.021 0.0239 1% Epsilon2c 14

31 0.014 0.020 0.015 0.022 0.030 0.030 0.0095 1];% Gamma12c 15

32

33 Q=data(var,mat);

34 end

39

C Funcao para o calculo da matriz de rigidez local e global q.m

1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


3 % Engenharia Mecânica− EESC USP



6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

7 function [Q]=q(i,layup);

8

9 mtrl=layup(i,1); %Reads the material number for the lamina i from the layup variable

10 teta=pi/180*layup(i,2); %Reads the angle for the lamina i from the layup variable

11

12 e1=data(2,mtrl); % Mechancal properties from database (data.m)

13 e2=data(3,mtrl);

14 v12=data(5,mtrl);

15 v21=e2*v12/e1;

16 g12=data(4,mtrl);

17

18 ql=1/(1−v12*v21)*[e1 e1*v21 0 % Local Stiffness Matrix

19 e2*v12 e2 0

20 0 0 g12*(1−v12*v21)];

21 c=cos(teta);

22 s=sin(teta);

23 T=[c^2 s^2 −2*c*s %Transformation matrix

24 s^2 c^2 2*c*s

25 c*s −c*s c^2−s^2];

26

27 Q=T*ql*T'; % Stiffness matrix in global coordinates

40

D Funcao para o calculo da matriz de rigidez local ql.m

1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

7 function [QL]=ql(i,layup);

8 mtrl=layup(i,1);

9 e1=data(2,mtrl);

10 e2=data(3,mtrl);

11 v12=data(5,mtrl);

12 v21=e2*v12/e1;

13 g12=data(4,mtrl);

14

15 QL=1/(1−v12*v21)*[e1 e1*v21 0

16 e2*v12 e2 0

17 0 0 g12*(1−v12*v21)];

18 end

41

E Funcao para o calculo da matriz de transformacao T.m

1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

7 function [t]=T(i,layup);

8 teta=pi/180*layup(i,2);

9 c=cos(teta);

10 s=sin(teta);

11 t=[c^2 s^2 −2*c*s

12 s^2 c^2 2*c*s

13 c*s −c*s c^2−s^2];

14 end

Documents

Programa Computacional Para Analise e Dimensionamento de Laminados Compositos Em MATLAB