Esercitazioni: travi reticolari, FEM calore 1-2D

TRASMISSIONE DEL CALORE MONODIMENSIONALE

RESIDUI PESATI - FUNZIONI DI FORMA

Nome: Federico

Cognome: Camiletti

Matricola: 0000805015

———————————————————————————————————————————

Introduzione 2

a) Residui Pesati 4

• Point collocation 6

• Subdomain collocation 9

Minimi quadrati 11

• Galerkin 13

•

b) Funzioni di forma 15

Lineari e paraboliche 18

Uso contemporaneo 26

1 di 26

Introduzione

Il metodo agli elementi finiti (FEM) è un metodo numerico che permette di ottenere una soluzione

approssimata di un generico problema differenziale esprimibile come:

(1) Lu(x) + g(x) = 0 dominio di definizione a<=x<=b e assegnate condizioni al contorno BC

L generico operatore differenziale, g(x) funzione nota e u(x) funzione incognita.

Il metodo si basa sui seguenti step:

Stabilire la formulazione forte del problema - Strong Form

• Ottenere la formulazione debole - Weak Form

• Scegliere un’approssimazione della funzione incognita - Uapp

• Scegliere la funzione di peso - v

•

Nell’esercitazione viene analizzata la trasmissione di calore monodimensionale in un elemento di

area A e conduttività termica k costanti, sottoposto a un carico noto Q(x). La (1) è definita come:

(2) STRONG FORM

In questo caso l’operatore L è d/dx^2, la funzione incognita T(x) e si hanno le due condizioni al

contorno (BC) T(a) e T(b). Moltiplicando la (2) per una funzione di peso arbitraria v, integrando sul

dominio e utilizzando l’integrazione per parti si ottiene la formulazione debole:

(3) WEAK FORM

La soluzione approssimata, scelta in modo da rispettare le BC, può essere espressa

come:

(4)

con: Parametri incogniti da determinare, in questo caso le temperature nei

nodi dell’elemento

Funzioni di forma che approssimano l’andamento della

temperatura nell’elemento

Il problema verrà risolto con con l’uso, anche contemporaneo, di diverse funzioni di forma.

2 di 26

Sostituendo in (1) si avrà

con errore e diverso da zero essendo la soluzione approssimata.

Moltiplicando per la funzione di peso v e integrando nel dominio si ha:

(5)

che rappresenta ancora la formulazione forte del problema

le funzioni di peso possono essere definite come

con: costanti arbitrarie funzioni di x note

In seguito verranno analizzati diversi metodi di residui pesati, che differiscono per la scelta

delle funzioni V.

Essendo c e a delle costanti si ha

e la (5) diventa:

Ponendo:

(6)

la risoluzione del sistema lineare Ka=f relativo alla formulazione debole del problema,

permette di ottenere ottenere le temperature nodali a e successivamente dalla (4) si

ottiene la l’andamento della soluzione approssimata. 3 di 26

a) Residui Pesati

Viene risolto il seguente problema: STRONG FORM (2)

Utilizzando per la scelta delle funzioni di peso V i metodi:

Point collocation

• Subdomain collocation

• Minimi quadrati

• Galerkin

•

Utilizzando Mupad si ottiene la soluzione esatta risolvendo separatamente due problemi

per le differenti condizioni di carico e imponendo la continuità nel punto intermedio. 4 di 26

Si può notare un andamento della funzione curvilineo e non simmetrico

Come soluzione approssimata si è scelto Uapp =

che rispetta le condizioni al contorno

in MuPAD

Funzioni di forma N B=L(N)

Dopo aver determinato la (6) per i diversi metodi si risolve il sistema lineare e viene

ricavata la soluzione approssimata. 5 di 26

• Point collocation

V si basa sulla funzione delta di Dirac e vale

In questo modo si impone che il residuo sia nullo in corrispondenza dei punti di

collocazione scelti e si ottiene:

delete K;

K := matrix(ind, ind):

for i from 1 to ind do

for j from 1 to ind do

K[i,j]:=float(subs(B[j,1],x=X[i,1]))

end_for

end_for:

K;

Si riportano la soluzione del problema ottenute per pochi punti di collocazione a diverse

posizioni:

per un punto di collocazione situato a metà intervallo in corrispondeva del quale si è posto

un carico pari a 0.5 si ottiene una soluzione simmetrica che non approssima bene quella

esatta. 6 di 26

Utilizzando due punti si ottiene un’approssimazione migliore, ma si può notare come

questa dipenda dalla posizione dei punti scelti; è importante quindi che i punti siano preso

il più possibile equidistanti e sufficientemente lontano dagli estremi.

x1 = 0.33 x2 = 0.66 x1 = 0.25 x2 = 0.75

Il codice è stato modificato per utilizzare un numero indefinito di punti, senza che vengano

presi in input le loro coordinate, prendendo il primo e l’ultimo punto a una distanza dagli

estremi pari a 1/20 della lunghezza e i punti intermedi equispaziati nell’intervallo:

Si può notare che selezionando solo 2 punti

vengono utilizzato gli estremi e per via della

disposizione non omogenea si ottiene una

soluzione per niente affidabile in particolare

lontano dai punti scelti. 7 di 26

Utilizzando già soli 7 punti si ottiene invece una buona approssimazione

Si può notare che scegliendo sempre 7 punti di collocazione, ma scegliendo anche gli

estremi x1=0 e xn=1 si ottiene un risultato peggiore 8 di 26

• Subdomain collocation

V viene definita come

In questo modo si impone che il residuo sia nullo in media in ogni sottodominio. Integrando

negli estremi degli intervalli si ottiene:

Dove nel codice di calcolo di F è stata considerata la diversa distribuzione del carico q,

nel sottodominio considerato. 9 di 26

Anche in questo caso per ottenere l’approssimazione migliore è opportuno scegliere

intervalli omogenei

[0,0.2] [0.2,1] [0,0.5] [0.5,1]

Il codice è stato modificato per utilizzare un numero indefinito di intervalli della stessa

lunghezza. Si riportano i risultati per 7 intervalli 10 di 26

Minimi quadrati

•

Il metodo impone che sia nullo il quadrato del residuo integrato sul tutto dominio.

V viene definita come

ed essendo

si ottiene:

con K simmetrica

Dove nel codice di calcolo di F è stata considerata la diversa distribuzione del carico q nell’

intero dominio. 11 di 26

In questo metodo è possibile utilizzare gli estremi dell’intervalli e si ottengono buoni

risultati anche con punti non equispaziati, inoltre si hanno vantaggi computazionali nella

risoluzione per via della simmetria di K. Si riportano i risultati per 7 punti equispaziati:

12 di 26

Galerkin

•

V viene definita come

essendo le funzioni di forma ortogonali al residuo

si ottiene: 13 di 26

Si riportano i risultati per 7 punti equispaziati

A parità di punti si ottiene la migliore approssimazione utilizzando questo metodo e inoltre, pur non

essendo in generale K simmetrica come nel metodo dei minimi quadrati, combinata alla

formulazione debole del problema permette di ottenere sempre una matrice simmetrica. Per questi

motivi il metodo di Galerkin è il più utilizzato per la scelta delle funzioni di peso nel metodo agli

elementi finiti. 14 di 26

b) Funzioni di forma

Viene risolto il problema utilizzando il metodo agli elementi finiti: STRONG FORM (2)

Si ottiene la seguente formulazione debole (3) dove, utilizzando il metodo di Galerking, le funzioni

di peso V vengono scelte uguali alle funzioni di forma N:

Dove in ordine troviamo: matrice di rigidezza K, vettore delle temperature nodali incognite a,

vettore delle condizioni al contorno e vettore dei carichi.

Nella risoluzione si utilizza il codice fornito FEMheat1d.m che, definiti i dati del problema, esegue i

seguenti passi:

discretizzazione dominio

• per ogni elemento:

• vettore dei carichi

• condizioni al contorno

• matrice di rigidezza

•

assemblaggio

• vincolamento

• temperature nodali a=K\f

• soluzione approssiomata

•

Il codice, che utilizza elementi con funzioni di forma lineari è stato modificato per utilizzare funzioni

paraboliche e di grado superiore e per permettere l’uso contemporaneo di diverse funzioni di

forma. 15 di 26

Per ricavare i diversi termini derivanti dalla funzione di forma N scelta si è creato un codice di

calcolo MuPAD phi_FEM.mn

dove per la funzione di forma N si utilizza la formula interpolatoria di Lagrange

per elementi di lunghezza L e definiti da nodi equispaziati in modo da ottenere un’espressione

semplificata dei termini da implementare in FEMheat1d.m. 16 di 26

LINEARE PARABOLICA

N

intN

K=N’N

Essendo la soluzione esatta del problema:

in seguito non saranno riportati i risultati per funzioni di forma di grado superiore.

Si riporta uno schema riassuntivo dei passaggi da eseguire nella risoluzione del sistema 17 di 26

Lineari e paraboliche

Dopo aver definito i dati del problema viene richiesto

in input il numero di elementi (di stessa lunghezza)

da generare e il loro grado.

Vengono calcolati i punti per elemento, creato il

vettore dei nodi equispaziali e la matrice di

connettività degli elementi. Utilizzando 3 elementi

lineari si ottengono: 18 di 26

Viene richiamata la funcion formStiffness1Dheat.m che tramite un ciclo su ogni elemento

determina il suo grado e calcola la sua matrice di rigidezza come K=ka/xa*C con conduttività

k ,area a, lunghezza xa dell’elemento e C matrice di coefficienti (ottenuta con MuPAD) scelta con

un comando switch in base al grado della funzione di forma dell’elemento.

Infine procede all’assemblaggio

con la matrice globale in

corrispondenza dei gradi di libertà

dell’elemento. 19 di 26

Procede allo stesso modo per il calcolo dei vettori dei carichi e delle condizioni al contorno:

In seguito al vincolamento

e dopo la correzione del vettore delle forze per via del vincolo imposto

la sottomatrice di K ottenuta eliminando riga e

colonna corrispondente al nodo vincolato

diventa invertibile e si ottiene la temperatura

nodale come soluzione del sistema tramite la

funcion: 20 di 26

Imponendo il valore nel nodo vincolato si ottiene:

mentre Ia funcion: permette di ricavare il flusso incognito al nodo 1 dovuto

all’eliminazione della della prima riga di K nel

vincolamento 21 di 26

l’andamento approssimato della temperatura sull’elemento viene tramite le relative funzioni di

forma e le temperature nodali: 22 di 26

Per ottenere una migliore approssimazione si utilizzano 7 elementi finiti,

sempre con funzioni di forma lineari, ottenendo i seguenti risultati: 23 di 26

Si riportano i risultati per 2 elementi finiti con funzioni di