Software – Modellazione magneti permanenti

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244

!     Last change:  PS   17 Apr 2003    5:48 pm
 
Autori: 
Paolo Sordi - Ordine degli ingegneri di Roma N° 22979 - Sezione A
Leonardo Santini - Ordine degli ingegneri di Roma N° 22757 - Sezione A
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE AZZERA_LA_MATRICE_GAUSSJ ()
USE CONFIG_FEM
USE CONFIG_MESH 
IMPLICIT NONE
INTEGER i,j
FUNZIONE_SCALARE = 0.0
VETTORE_PRODOTTO = 0.0
DO i = 1, NNOD
    DO j = 1,NNOD
        MATRICE_GLOBALE_S (i,j) = 0.0
    END DO
END DO
END SUBROUTINE
 
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE IMPOSTA_CONDIZIONI_INIZIALI_NEL_SISTEMA_LINEARE_GAUSSJ ()
USE WINTERACTER
USE CONFIG_FEM
USE CONFIG_MESH 
IMPLICIT NONE
INTEGER :: k
INTEGER :: i
DO k = 1 , NNOD
   IF (CONDIZIONE_INIZIALE(k).NE.-1.0) THEN
       DO i=1,NNOD
           MATRICE_GLOBALE_S (k,i) = 0.0
       END DO
       MATRICE_GLOBALE_S (k,k) = 1.0
       VETTORE_PRODOTTO (k) = CONDIZIONE_INIZIALE (k)
       WRITE (1000,*) 'CONDIZIONE_INIZIALE (k)', CONDIZIONE_INIZIALE (k)
   END IF
END DO
PAUSE
END SUBROUTINE
 
 
 
!*************************************************************************
!*************************************************************************
 
 
 
SUBROUTINE CREAZIONE_MATRICE_S_PER_MATERIALI_LINEARI_gaussj ()
USE CONFIG_MESH
USE CONFIG_FEM
USE WINTERACTER
IMPLICIT NONE
DOUBLE PRECISION,DIMENSION (3)   :: grad1,grad2,grad3,grad4
DOUBLE PRECISION,DIMENSION (4,4) :: Matrix_locale_S
DOUBLE PRECISION                 :: dot
DOUBLE PRECISION                 :: value
DOUBLE PRECISION                 :: CALCOLO_VOLUME_TETRAEDRO
DOUBLE PRECISION                 :: permeabilita_relativa
real, DIMENSION (4,4)            :: matrix_ABCD
INTEGER i,k,k1
INTEGER irow,icol
 
DO k= 1,NTETRA
    IF (MATERIALE_TETRAEDRO (k).eq.1.or.MATERIALE_TETRAEDRO (k).eq.2) THEN
        CALL CALCOLA_COEFFICIENTI_A_B_C_D (k, matrix_ABCD)
        DO i = 1,3
            grad1 (i) = matrix_ABCD (i+1,1)
            grad2 (i) = matrix_ABCD (i+1,2)
            grad3 (i) = matrix_ABCD (i+1,3)
            grad4 (i) = matrix_ABCD (i+1,4)
        END DO
        Matrix_locale_S = 0.0
        Matrix_locale_S (1,1) = DOT (grad1,grad1)
        Matrix_locale_S (2,2) = DOT (grad2,grad2)
        Matrix_locale_S (3,3) = DOT (grad3,grad3)
        Matrix_locale_S (4,4) = DOT (grad4,grad4)
        Matrix_locale_S (1,2) = DOT (grad1,grad2)
        Matrix_locale_S (1,3) = DOT (grad1,grad3)
        Matrix_locale_S (1,4) = DOT (grad1,grad4)
        Matrix_locale_S (2,3) = DOT (grad2,grad3)
        Matrix_locale_S (2,4) = DOT (grad2,grad4)
        Matrix_locale_S (3,4) = DOT (grad3,grad4)
!       CALL CONTROLLA_CHE_LA_MATRICE_SIA_SINGOLARE (Matrix_locale_S)
!        Matrix_locale_S = Matrix_locale_S  / (36.0 * CALCOLO_VOLUME_TETRAEDRO(k))
        Matrix_locale_S = Matrix_locale_S * CALCOLO_VOLUME_TETRAEDRO(k)
        IF (MATERIALE_TETRAEDRO (k).eq.1) THEN
            permeabilita_relativa = 1.0
            Matrix_locale_S = Matrix_locale_S * (2.  * permeabilita_relativa)
        END IF
        IF (MATERIALE_TETRAEDRO (k).eq.2) THEN
            permeabilita_relativa = 100
            Matrix_locale_S = Matrix_locale_S * (2.  * permeabilita_relativa)
        END IF
        DO k1= 1,4
            MATRICE_GLOBALE_S (NT(k,k1),NT(k,k1)) = MATRICE_GLOBALE_S (NT(k,k1),NT(k,k1)) + Matrix_locale_S (k1,k1)
        END DO
        MATRICE_GLOBALE_S (NT(k,1),NT(k,2)) = MATRICE_GLOBALE_S (NT(k,1),NT(k,2)) + Matrix_locale_S (1,2)
        MATRICE_GLOBALE_S (NT(k,1),NT(k,3)) = MATRICE_GLOBALE_S (NT(k,1),NT(k,3)) + Matrix_locale_S (1,3)
        MATRICE_GLOBALE_S (NT(k,1),NT(k,4)) = MATRICE_GLOBALE_S (NT(k,1),NT(k,4)) + Matrix_locale_S (1,4)
        MATRICE_GLOBALE_S (NT(k,2),NT(k,3)) = MATRICE_GLOBALE_S (NT(k,2),NT(k,3)) + Matrix_locale_S (2,3)
        MATRICE_GLOBALE_S (NT(k,2),NT(k,4)) = MATRICE_GLOBALE_S (NT(k,2),NT(k,4)) + Matrix_locale_S (2,4)
        MATRICE_GLOBALE_S (NT(k,3),NT(k,4)) = MATRICE_GLOBALE_S (NT(k,3),NT(k,4)) + Matrix_locale_S (3,4)
 
        MATRICE_GLOBALE_S (NT(k,2),NT(k,1)) = MATRICE_GLOBALE_S (NT(k,2),NT(k,1)) + Matrix_locale_S (1,2)
        MATRICE_GLOBALE_S (NT(k,3),NT(k,1)) = MATRICE_GLOBALE_S (NT(k,3),NT(k,1)) + Matrix_locale_S (1,3)
        MATRICE_GLOBALE_S (NT(k,4),NT(k,1)) = MATRICE_GLOBALE_S (NT(k,4),NT(k,1)) + Matrix_locale_S (1,4)
        MATRICE_GLOBALE_S (NT(k,3),NT(k,2)) = MATRICE_GLOBALE_S (NT(k,3),NT(k,2)) + Matrix_locale_S (2,3)
        MATRICE_GLOBALE_S (NT(k,4),NT(k,2)) = MATRICE_GLOBALE_S (NT(k,4),NT(k,2)) + Matrix_locale_S (2,4)
        MATRICE_GLOBALE_S (NT(k,4),NT(k,3)) = MATRICE_GLOBALE_S (NT(k,4),NT(k,3)) + Matrix_locale_S (3,4)
   END IF
END DO
CALL WMESSAGEBOX (0,0,1,'Variabili azzerate. Fine processo di immissione della matrice','')
END SUBROUTINE
 
 
 
!*************************************************************************
!*************************************************************************
 
 
 
SUBROUTINE CREAZIONE_MATRICE_S_E_VETTORE_PER_MAGNETE_PERMANENTE_GAUSSJ ()
USE CONFIG_MESH
USE CONFIG_FEM
USE WINTERACTER
IMPLICIT NONE
DOUBLE PRECISION,DIMENSION (3)   :: grad1,grad2,grad3,grad4
DOUBLE PRECISION,DIMENSION (4,4) :: Matrix_locale_S
DOUBLE PRECISION                 :: dot
DOUBLE PRECISION                 :: value
DOUBLE PRECISION :: CALCOLO_VOLUME_TETRAEDRO
DOUBLE PRECISION :: Campo_coercitivo_Hc
DOUBLE PRECISION :: Br_componente_x
DOUBLE PRECISION :: Br_componente_y
DOUBLE PRECISION :: Br_componente_z
DOUBLE PRECISION :: Br_induzione_residua
DOUBLE PRECISION :: componente_vettore
real, DIMENSION (4,4) :: matrix_ABCD
INTEGER i,k,k1
INTEGER irow,icol
! CALCOLO DEGLI ELEMENTI DELLA MATRICE APPARTENENTI AL MAGNETE
OPEN  (4000,FILE= "componente magnete matrice_superiore.txt")
 
Campo_coercitivo_Hc  = 800000. ! Ampere/m
Br_componente_x = 1.2
Br_componente_y = 0.0
Br_componente_z = 0.0
Br_induzione_residua = DSQRT (Br_componente_x ** 2.0 + Br_componente_y ** 2.0 + Br_componente_z ** 2.0)
DO k = 1, NTETRA
   IF (MATERIALE_TETRAEDRO (k).EQ.3) THEN
       CALL CALCOLA_COEFFICIENTI_A_B_C_D (k, matrix_ABCD)
       DO i = 1,3
           grad1 (i) = matrix_ABCD (i+1,1)
           grad2 (i) = matrix_ABCD (i+1,2)
           grad3 (i) = matrix_ABCD (i+1,3)
           grad4 (i) = matrix_ABCD (i+1,4)
       END DO
       Matrix_locale_S = 0.0
       Matrix_locale_S (1,1) =  DOT(grad1,grad1)
       Matrix_locale_S (2,2) =  DOT(grad2,grad2)
       Matrix_locale_S (3,3) =  DOT(grad3,grad3)
       Matrix_locale_S (4,4) =  DOT(grad4,grad4)
       Matrix_locale_S (1,2) =  DOT(grad1,grad2)
       Matrix_locale_S (1,3) =  DOT(grad1,grad3)
       Matrix_locale_S (1,4) =  DOT(grad1,grad4)
       Matrix_locale_S (2,3) =  DOT(grad2,grad3)
       Matrix_locale_S (2,4) =  DOT(grad2,grad4)
       Matrix_locale_S (3,4) =  DOT(grad3,grad4)
!       CALL CONTROLLA_CHE_LA_MATRICE_SIA_SINGOLARE (Matrix_locale_S)
       Matrix_locale_S = Matrix_locale_S * (Campo_coercitivo_Hc / Br_induzione_residua)
!       Matrix_locale_S = Matrix_locale_S / (36.0 * CALCOLO_VOLUME_TETRAEDRO(K))
       Matrix_locale_S = Matrix_locale_S * CALCOLO_VOLUME_TETRAEDRO(K)
       Matrix_locale_S = Matrix_locale_S / 2.0
       DO k1= 1,4
            MATRICE_GLOBALE_S (NT(k,k1),NT(k,k1)) = MATRICE_GLOBALE_S (NT(k,k1),NT(k,k1)) + Matrix_locale_S (k1,k1)
       END DO
       MATRICE_GLOBALE_S (NT(k,1),NT(k,2)) = MATRICE_GLOBALE_S (NT(k,1),NT(k,2)) + Matrix_locale_S (1,2)
       MATRICE_GLOBALE_S (NT(k,1),NT(k,3)) = MATRICE_GLOBALE_S (NT(k,1),NT(k,3)) + Matrix_locale_S (1,3)
       MATRICE_GLOBALE_S (NT(k,1),NT(k,4)) = MATRICE_GLOBALE_S (NT(k,1),NT(k,4)) + Matrix_locale_S (1,4)
       MATRICE_GLOBALE_S (NT(k,2),NT(k,3)) = MATRICE_GLOBALE_S (NT(k,2),NT(k,3)) + Matrix_locale_S (2,3)
       MATRICE_GLOBALE_S (NT(k,2),NT(k,4)) = MATRICE_GLOBALE_S (NT(k,2),NT(k,4)) + Matrix_locale_S (2,4)
       MATRICE_GLOBALE_S (NT(k,3),NT(k,4)) = MATRICE_GLOBALE_S (NT(k,3),NT(k,4)) + Matrix_locale_S (3,4)
 
       MATRICE_GLOBALE_S (NT(k,2),NT(k,1)) = MATRICE_GLOBALE_S (NT(k,2),NT(k,1)) + Matrix_locale_S (1,2)
       MATRICE_GLOBALE_S (NT(k,3),NT(k,1)) = MATRICE_GLOBALE_S (NT(k,3),NT(k,1)) + Matrix_locale_S (1,3)
       MATRICE_GLOBALE_S (NT(k,4),NT(k,1)) = MATRICE_GLOBALE_S (NT(k,4),NT(k,1)) + Matrix_locale_S (1,4)
       MATRICE_GLOBALE_S (NT(k,3),NT(k,2)) = MATRICE_GLOBALE_S (NT(k,3),NT(k,2)) + Matrix_locale_S (2,3)
       MATRICE_GLOBALE_S (NT(k,4),NT(k,2)) = MATRICE_GLOBALE_S (NT(k,4),NT(k,2)) + Matrix_locale_S (2,4)
       MATRICE_GLOBALE_S (NT(k,4),NT(k,3)) = MATRICE_GLOBALE_S (NT(k,4),NT(k,3)) + Matrix_locale_S (3,4)
   END IF
END DO
CALL WMESSAGEBOX (0,0,1,'Fine del processo di immissione nella matrice S ','ciao')
 
DO k = 1 , NTETRA
   IF (MATERIALE_TETRAEDRO (k).EQ.3) THEN
       CALL CALCOLA_COEFFICIENTI_A_B_C_D (k, matrix_ABCD)
       DO i=1,4
             componente_vettore = Br_componente_x * matrix_ABCD (2,i)  &
                                + Br_componente_y * matrix_ABCD (3,i)  &
                                + Br_componente_z * matrix_ABCD (4,i)
             componente_vettore = ( Campo_coercitivo_Hc / Br_induzione_residua ) * componente_vettore
!            componente_vettore = componente_vettore / 2.0
             componente_vettore = componente_vettore * CALCOLO_VOLUME_TETRAEDRO (k)
             VETTORE_PRODOTTO (NT(k,i)) = VETTORE_PRODOTTO (NT(k,i)) + componente_vettore
       END DO
   END IF
END DO
CALL WMESSAGEBOX (0,0,1,'Fine del processo di immissione nel vettore forzamento','')
END SUBROUTINE
 
 
!*************************************************************************
!*************************************************************************
 
 
SUBROUTINE SOLUZIONE_ELEMENTI_FINITI_CON_GAUSSJ
USE CONFIG_FEM
USE CONFIG_MESH
REAL ,DIMENSION (MAXNODI) :: VETT_PRODOTTO
WRITE (1000,*) "MATRICE_GLOBALE_S"
DO i = 1, NNOD
    DO j = 1,NNOD
        IF (MATRICE_GLOBALE_S (i,j).ne.0.0) WRITE (1000,*) i,j, MATRICE_GLOBALE_S (i,j)
    END DO
END DO
WRITE (1000,*)  'VETT_PRODOTTO'
VETT_PRODOTTO = REAL (VETTORE_PRODOTTO)
DO j = 1,NNOD
    WRITE (1000,*) j, VETT_PRODOTTO (j)
END DO
CALL GAUSSJ (MATRICE_GLOBALE_S,NNOD,MAXNODI,VETT_PRODOTTO,1,1)
FUNZIONE_SCALARE = VETTORE_PRODOTTO
PAUSE 
END SUBROUTINE
On: